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Abstract . 

We  have  used  data  from  the  Mark  III  very-long-baseline  interferometry  (VLBI) 
system  to  study  the  relative  motions  of  radiotelescopes  on  the  Pacific,  North  Ameri¬ 
can,  and  Eurasian  plates,  and  in  the  deformation  zone  between  the  North  American 
and  Pacific  plates  in  the  California  region.  These  results  are  in  accord  with  recent 
geologic  plate  motion  models,  and  a  distributed  deformation  zone  between  the  North 
American  and  Pacific  plates.  We  have  carried  out  a  number  of  studies  on  the  accu¬ 
racy  of  the  results  being  obtained  with  VLBI.  The  most  notable  of  these  studies  are 
of  the  repeatability  of  the  estimates  of  three-dimensional  coordinates  of  the  sites  for 
long  baselines  (&4,000  km),  and  the  effects  of  calibrations  of  the  “wet”  tropospheric 
delays  using  water-vapor  radiometers.  We  have  also  been  studying  the  usefulness  of 
observations  made  at  very  low  elevation  angles  (r$4/*)  for  improving  the  accuracy  of 
VLBI  baseline  determinations.  We  have  also  studied  the  forced  nutations  of  the  Earth 
which  yielded  results  consistent  with  the  flattening  of  the  core-mantle  boundary  being 
about  5/o  greater  than  that  predicted  by  the  assumption  of  hydrostatic  equilibrium 
within  the  Earth.  The  effects  of  ocean  tides  and  the  anelasticity  of  the  mantle  on  the 
forced  nutations  are  not  clearly  evident  in  our  results,  and  are  still  being  investigated. 
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Significant  Accomplishments 

Under  this  contract  we  have  been  examining  many  facets  of  the  geodetic  and 
geophysical  applications  of  very-long-baseline  interferometry  (VLBI).  In  the  following 
sections  we  present  a  summary  of  the  results  which  have  been  obtained.  In  the  ap¬ 
pendices  to  this  report  we  present  a  detailed  discussion  of  the  procedures  we  are  using 
in  our  Kalman  filter  VLBI  data  analysis  software,  and  we  include  pre-  and  reprints  of 
the  papers  we  have  published  under  this  contract  (Davis  et  al .,  1988;  Elgered  et  al., 
1988;  Herring  et  al .,  1988;  Herring  ,1988) 


j 

I 


Global  site  velocity  estimates 


For  our  analysis  of  the  VLBI  data  set,  we  combined  532  experiments  carried 
out  between  July  1980  and  June  1987.  This  data  set  contained  202,051  group  delay 
measurements.  From  21,245  estimates  of  site  positions,  189  global  parameters  were 
estimated.  We  restricted  the  sites  used  in  this  analysis  to  those  which  had  been  used 
five  or  more  times  during  the  last  seven  years.  The  global  parameters  consisted  of  31 
site  position  estimates  and  30  site  velocity  estimates.  (One  site  was  used  for  only  a 
one  week  period,  and,  hence,  no  velocity  was  estimated  for  this  site).  We  estimated 
all  site  positions  and  velocities  (except  as  noted  above)  with  weak  a  priori  constraints 
applied  to  these  quantities  (±1  meter  for  positions,  and  ±1  meter/year  for  velocities). 
In  addition,  we  allowed  the  position  of  the  earth’s  rotation  axis  and  the  value  of 
UTl-AT  to  change  randomly  between  experiments.  When  all  of  the  experiments  are 
combined  in  this  fashion,  the  final  estimates  of  the  site  position  and  velocities  have 
large  absolute  uncertainties  because  the  coordinate  system  in  which  they  are  given  is 
not  well  defined.  Since  all  sites  were  allowed  to  move,  and  the  changes  in  the  earth 
orientation  parameters  were  weakly  constrained,  the  coordinate  system  in  which  the 
final  results  are  given  can  translate  and  rotate  relative  to  its  position  and  orientation 
at  the  epoch  of  the  first  experiment  in  the  solution.  To  remove  this  rank  deficiency  in 
the  final  results  we  need  to  adopt  the  values  for  three  scalars  to  represent  a  general 
translation,  and  three  scalars  to  represent  a  general  rotation  of  the  coordinate  system, 
relative  to  its  initial  position  and  orientation.  There  are  no  unique  values  for  these 
six  parameters,  although  there  are  a  number  of  methods  we  could  adopt  for  obtaining 
values  of  these  parameters.  For  example,  the  motion  of  one  site  could  be  adopted  as  a 
reference  in  order  to  determine  the  translation  of  the  coordinate  system,  and  the  earth 
orientation  values  from  an  earth  orientation  service  (such  a*  BIH  or  IRIS)  could  be 
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adopted  to  define  the  final  orientation  of  the  coordinate  system.  However,  these  types 
of  techniques  do  not  provide  a  very  satisfactory  solution  to  the  problem.  If  the  motion 
of  the  reference  is  different  from  the  adopted  motion,  then  the  difference  will  be  applied  I 

to  all  other  sites  in  the  network,  thus  leading  to  apparent  motions  at  these  sites.  The 
adoption  of  earth  orientation  values  from  BIH  or  IRIS  also  presents  problems  because 
neither  of  these  services  currently  includes  site  motions  in  its  determinations  of  the 
earth  orientation  parameters.  ^ 

The  procedure  we  have  adopted  for  determining  the  three  translation  and  three 
rotation  parameters  is  to  estimate  their  values  by  minimizing  the  root-mean-square 
(RMS)  difference  between  the  estimated  (locally)  horizontal  velocities  and  the  pre- 
dieted  velocities  from  a  global  plate  motion  model  for  a  selected  set  of  sites.  To  resolve 
the  rank  deficiency  for  this  initial  solution,  we  used  the  Minster  and  Jordan  RM2  plate 
motion  model,  and  the  velocities  of  11  of  the  30  sites  whose  velocities  were  estimated. 

These  eleven  sites  were  selected  because  of  their  geographical  coverage  of  the  VLBI  I 

network,  and  the  frequency  of  their  use  during  the  seven  years  spanned  by  the  VLBI 
data  set.  The  estimated  velocties  and  their  one  sigma  standard  deviations,  and  the 
RM2  velocities  for  all  30  sites,  are  given  in  Table  1.  These  results  are  also  presented 
in  graphical  form  in  Figures  1-3.  The  figures  clearly  show  a  pattern  of  motion  which  j 

is  very  consistent  with  the  RM2  plate  motion  model  except  in  the  transition  region 
between  the  North  American  and  Pacific  plates  along  the  San  Andreas  fault.  Even 
here,  the  motions  seem  to  show  a  smooth  change  between  North  American  velocities 
and  Pacific  velocities.  | 

There  were  some  unexpected  results  which  came  out  of  this  solution.  The  relative 
(horizontal)  motions  of  the  VLBI  sites  in  Eastern  North  America,  including  Ft.  Davis, 

TX,  are  all  very  small,  and  consistent  with  zero  within  their  90%  confidence  intervals. 

This  result  was  surprising  because  the  length  of  the  baseline  between  the  Westford,  ! 

MA,  and  the  Ft.  Davis,  TX,  sites  has  been  decreasing  on  average  by  ss7  mm/yr  for  the 

last  seven  years,  although  there  have  been  a  number  of  excursions  from  this  average 

rate.  For  the  last  year,  the  length  has  shown  no  appreciable  change,  and  during  this 

interval  the  length  is  very  similar  to  its  value  in  1980.  In  our  global  solution,  this  1 

decreasing  baseline  length  is  due  mainly  to  vertical  motions  at  the  Ft.  Davis  site.  (For 

the  Westford-Ft.  Davis  baseline,  «25%  of  the  height  change  projects  into  the  baseline 

length.)  We  are  still  investigating  whether  this  height  change  is  due  to  the  motion  of 

the  site,  or  due  to  an  inadequacy  in  our  models  of  the  tropospheric  delay.  All  of  the  ' 
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TABLE  1.  Estimated  velocities  of  VLBI  sites  obtained  from  the  analysis  discussed  in 
the  text. 


Site 

<t> 

(deg) 

A 

(deg) 

Estimates  RM2  Predicted 

North  Vel.  East  Vel.  North  Vel.  East  Vel. 
(mm/yr)  (mm/yr)  (mm/yr)  (mm/yr) 

Plate* 

Haystack*,MA 

42.6 

288.5 

-5. ±2. 

-22. ±3. 

-4. 

-24. 

NA 

Westford*,MA 

42.6 

288.5 

-5. ±2. 

-23. ±3. 

-4. 

-24. 

NA 

Maryland  Pt.,MD 

38.4 

282.8 

-8. ±9. 

-19. ±8. 

-5. 

-24. 

NA 

NRAO,WV 

38.4 

280.2 

-13. ±6. 

-28. ±5. 

-6 

-24. 

NA 

Richmond, FL 

25.6 

279.6 

-10. ±3. 

-25. ±3. 

-6. 

-25. 

NA 

Ft.  Davis, TX 

30.6 

256.1 

-8. ±2. 

-23. ±2. 

-9. 

-23. 

NA 

Platville,CO 

40.2 

255.3 

-19. ±6. 

-19. ±5. 

-9. 

-22. 

NA 

Flagstaff,  A  Z 

35.2 

248.4 

-6. ±10. 

-39. ±9. 

-10. 

-21. 

NA 

Yuma,AZ 

32.9 

245.8 

-10. ±3. 

-29. ±3. 

-10. 

-21. 

NA 

Ely,NV 

39.3 

245.6 

-16. ±5. 

-27. ±4. 

-10. 

-20. 

NA 

Blackbutte.CA 

33.7 

244.3 

19. ±6. 

-42. ±5. 

-10. 

-21. 

NA 

Monument  Peak, C A 

32.9 

243.6 

19. ±3. 

-50. ±2. 

-10. 

-21. 

NA 

Pinyon  Flats, CA 

33.6 

243.5 

17. ±4. 

-46. ±3. 

-10. 

-21. 

NA 

Mojave*,  CA 

35.3 

243.1 

-3. ±2. 

-27. ±2. 

-10. 

-21. 

NA 

Pearblossom,CA 

34.5 

242.1 

8. ±7. 

-38. ±5. 

-10. 

-21. 

NA 

Pasadena, CA 

34.2 

241.8 

l.±3. 

-46. ±2. 

-10. 

-21. 

NA 

Owens  Valley*,CA 

37.2 

241.7 

-4. ±2. 

-28. ±2. 

-10. 

-21. 

NA 

Mammouth  Lakes, CA 

37.6 

241.1 

-3. ±5. 

-30. ±3. 

-10. 

-21. 

NA 

Vandenberg*,CA 

34.6 

239.3 

26. ±2. 

-56. ±2. 

35. 

-55. 

PA 

Quincy, CA 

40.0 

239.1 

-4. ±3. 

-28. ±3. 

-10. 

-20. 

NA 

Hatcreek*,CA 

40.8 

238.5 

-3. ±2. 

-26. ±2. 

-10. 

-20. 

NA 

Ft.  Ord.CA 

36.7 

238.2 

28. ±6. 

-49. ±5. 

-10. 

-21. 

NA 

Presidio, CA 

37.8 

237.5 

15. ±8. 

-42. ±6. 

-10. 

-21. 

NA 

Gilmore  Ck.*,AK 

65.0 

212.5 

-15. ±3. 

-8. ±3. 

-10. 

-8. 

NA 

Kauai*, HI 

22.1 

200.3 

44. ±5. 

-79.-4. 

53. 

-85. 

PA 

Kwajalean* 

9.4 

167.5 

53. ±9. 

-90.  ±8. 

49. 

-97. 

PA 

Kashima*,  Japan 

36.0 

140.1 

9. ±8. 

-21. ±8. 

O 

•  ar . 

-13. 

NA 

Wettzell*,  FRG 

49.1 

12.9 

2. ±3. 

5. ±4. 

0. 

0. 

EU 

Onsala*,  Sweden 

57.4 

11.9 

0.±3. 

-l.±4. 

0. 

0. 

EU 

*  The  abbreviations  of  the  tectonic  plate  names  are  NA-North  America,  PA-Pacific,  and 
EU-Eurasia.  The  velocities  are  arbitrarily  referred  to  the  Eurasian  plate. 

*  These  sites  were  used  in  the  determination  of  the  translation  and  the  orientation  of  the 
coordinate  system  (see  text  for  discussion) . 
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Figure  1.  Vector  motions  of  the  VLB!  sites  given  in  Table  I.  The  error  ellipses  give  90%  confidence  limits. 
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r  motions  of  the  VLHI  sites  given  in  Table  1  with  the  RM2  plate  motion  model  removed.  (Only  selected  sites 
e  region  of  California,  see  figure  .l.).  The  error  ellipses  give  90%  confidence  limits.  Note  the  change  in  scale 


other  frequently  observed  VLBI  sites  have  vertical  velocities  of  less  than  10  mm/yr 
with  standard  deviations  of  typically  a  few  mm/yr.  Our  analysis  also  shows  that  the 
VLBI  sites  at  Mojave,  Owens  Valley,  Hatcreek,  Mammouth  Lakes,  and  Quincy  are 
all  moving  relative  to  North  America  with  velocities  which  are  consistent  with  the 
extension  of  the  Basin  and  Range  Province  of  California. 

Three-dimensional  coordinates 

In  the  past  we  have  not  used  three-dimensional  coordinates  for  tectonic-motion 
studies  because  the  coordinate  system  in  which  these  values  are  determined  is  fixed 
with  respect  to  the  rotation  axis  of  the  Earth.  This  axis  moves  relative  to  the  crust  of 
the  Earth.  The  motion  of  the  rotation  axis  can  be  as  large  as  12  meters  during  a  12 
month  interval,  and  hence  must  be  accurately  accounted  for  in  any  study  of  changes 
in  the  three-dimensional  coordinates  of  VLBI  sites.  Conventional  geodetic  techniques 
are  not  accurate  enough  for  us  to  be  able  to  use  their  measurements  of  the  position 
of  the  Earth’s  rotation  axis  for  tectonic-m ution  studies.  However,  now  that  modern 
geodetic  techniques,  such  as  very-long-baseline  interferometry  (VLBI),  satellite  laser 
ranging  (SLR)  and  lunar  laser  ranging  (LLR),  are  being  used  to  monitor  the  motions 
of  the  Earth’s  rotation  axis,  we  can  now  consider  using  three-dimensional  coordinates 
for  tectonic-motion  studies. 

To  test  the  accuracy  of  the  Earth  rotation  measurements  obtained  from  modern 
techniques  we  have  been  studying  the  repeatability  of  the  estimates  of  the  coordinates 
of  the  VLBI  site  in  the  Mojave  desert  obtained  from  the  VLBI  experiments  being  used 
to  test  the  usefulness  of  observations  taken  at  low-elevations  angles.  These  experi¬ 
ments  are  often  referred  to  as  the  “LOWEL”  experiments.  We  have  concentrated  on 
comparing  the  results  from  two  different  sets  of  Earth  rotation  measurements.  One 
series  of  measurements,  which  we  will  refer  to  as  the  “CfA”  series,  was  obtained  from 
our  analysis  of  all  VLBI  data  except  for  those  in  the  LOWEL  experiments.  The  CfA 
series  is  a  North  America  fixed  system  (t  e.,  in  this  system,  sites  not  moving  relative 
to  the  North  American  plate,  should  have  constant  coordinates).  The  other  series, 
which  we  will  refer  to  as  the  ”JPL  PRIME”  series,  was  obtained  at  the  Jet  Propulsion 
Laboratory  by  combining  Earth  rotation  data  from  Satellite  Laser  Ranging,  Lunar 
Laser  Ranging  and  Deep-Space  Network  VLBI  data.  The  JPL  PRIME  series  should 
have  no  rotation  relative  to  the  AMO-1  global  plate  model.  Within  this  system,  sites 
on  the  North  American  plate  should  move  with  a  velocity  of  about  2  cm/yr.  These  two 
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series  of  Earth  rotation  values  contain  no  common  data.  In  Figures  4  and  5  we  show 
the  estimated  positions  of  the  Mojave  site  in  the  coordinate  systems  defined  by  these 
earth  rotation  series.  The  coordinates  are  shown  in  a  (Mojave)  local  north,  east  and  up 
(height)  frame.  The  repeatability  of  the  estima  .es  A  the  Mojave  coordinates  is  better 
with  the  CfA  earth  rotation  values  than  with  the  JPL  Prime  series.  The  vdocity  of  the 
Mojave  site  obtained  when  the  CfA  earth  rotation  series  is  used  is  consistent,  within 
the  uncertainties  of  these  estimates,  with  the  estimates  of  the  Mojave  velocity  obtained 
from  the  global  solutions  discussed  in  previous  reports.  (These  global  solutions  do  not 
include  the  LOWEL  experiments.  The  velocity  estimates  from  these  global  solutions 
range  between  0.9-0. 7  N  and  -0.7  and  -0.4  E  cm/yr,  depending  on  which  sites  are  used 
to  obtain  the  reference  frame;  see  previous  reports.) 

We  may  use  these  results  to  assess  the  precision  of  the  CfA  and  JPL  PRIME  earth 
rotation  series.  Most  of  the  LOWEL  experiments  have  been  single  baseline  experiments 
using  the  Haystack  Observatory  or  the  Westford  Observatory  as  the  second  site  in  the 
baseline.  Since  Haystack  and  Westford  are  only  1.24  kms  apart,  the  sensitivity  of 
the  Mojave  coordinates  to  changes  in  the  Earth  rotation  parameters  will  be  almost 
identical  for  these  two  baselines.  When  coordinates  of  either  Haystack  or  Westford  are 
fixed  at  their  apriori  values  (to  define  the  translation  origin  of  the  coordinate  system), 
the  sensitivity  of  the  Mojave  North,  East  and  height  estimates  ( 6N ,  6E  and  6U)  to 
changes  in  the  Earth  rotation  parameters,  6x,  6 y,  and  6UT\,  is  given  by: 

(6N\  /  — 14  -1  9  W  6x  \ 

6E  =  3  1  -9  6y  (1) 

\SU  J  \  —12  3  -13 )  \6UTl 

•where  the  coefficients  for  the  matrix  are  computed  for  coordinate  changes  in  mm  and 
the  orientation  changes  in  mas.  In  addition,  any  errors  in  modeling  the  atmospheric 
delays  at  the  Haystack/ Westford  site,  which  would  normally  affect  the  estimate  of  the 
height  at  these  sites,  will  project  into  the  coordinates  of  the  Mojave  site.  For  a  height 
error  of  6Ujj/w  at  these  sites,  the  effect  on  the  estimates  of  the  Mojave  coordinates  is 
be  given  by: 

/  0.25  \ 

SE  —  0.53  6Uf{/w.  (2) 

\6U  J  90  I 
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Figure  4.  Estimated  positions  of  the  Mojave  antenna  from  the  LOWEL  experiments 
•with  the  orientation  of  the  coordinate  system  defined  by  the  CfA  Earth  rotation  series 
a  (see  text).  The  coordinate  system  is  a  (Mojave)  North,  East  and  Up  frame,  where  N 

is  the  ellipsoidal  arc  distance  from  the  equator,  E  is  the  distance  from  the  Greenwich 
meridian  along  the  small  circle  at  the  geodetic  latitude  of  the  site,  and  U  is  the  height 
of  the  site  above  the  ellipsoid.  The  error  bars  do  not  account  for  the  uncertainty  in 
■  the  realization  of  the  North  American  fixed  coordinate  system. 
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Figure  5.  Estimated  positions  of  the  Mojave  antenna  from  the  LOWEL  experiments 
with  the  orientation  of  the  coordinate  system  defined  by  the  JPL  PRIME  Earth  rota¬ 
tion  series  (see  text).  See  Figure  1  caption  for  definition  of  N,  E  and  U. 
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If  we  assume  that  errors  in  the  estimates  of  the  North  and  East  components  of  the 
Mojave  site  position  are  dominated  by  errors  in  the  Earth  rotation  data  and  by  errors 
in  the  height  of  Haystack/VVestford,  then  we  may  use  Equations  1  and  2  to  estimate 
uncertainties  of  the  Earth  orientation  values  and  6Uh/w  Such  an  analysis  of  the  CfA 
results  yields,  assuming  that  the  variances  of  the  errors  in  each  of  the  components  of 
the  Earth  rotation  parameters  are  equal,  uncertainties  of  0.4  mas  for  the  Earth  rotation 
parameters  and  28  mm  for  6Uh/w.  A  similar  analysis  for  the  results  obtained  using 
the  JPL  PRIME  series,  yields  1.2  mas  for  the  Earth  rotation  parameters  and  40  mm 
for  6Uh/w ■  This  uncertainty  of  the  JPL  PRIME  series  is  consistent  with  the  quoted 
uncertainties  for  this  series,  and  is  confirmed  by  the  direct  comparison  of  the  CfA  and 
JPL  PRIME  series  shown  in  Figure  6. 

We  plan  to  study  further  the  usefulness  of  three-dimensional  coordinates  for 
studying  tectonic  motions.  We  will  probably  next  process  the  Onsala-Wettzell  (sslOOO 
km)  baseline  using  the  CfA  and  the  JPL  PRIME  earth  rotation  series,  and  assesses 
the  precision  with  which  differential  position  of  these  two  sites  can  be  determined. 


II 
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Appendix  A.  Kalman  filter  estimator 

We  discuss  here  the  basic  set  of  matrix  operations  needed  to  implement  a  Kalman 
filter  estimator.  These  equations  appear  in  most  text  book  which  discuss  Kalman 
filtering.  The  particular  form  we  use  is  based  on  the  works  of  Liebelt  (1967)  and  Gelb 
(1974).  We  commence  with  the  linearized  form  of  the  equations  which  relate  VLBI 
observables  to  the  parameters  to  be  estimated: 

yt  —  Atit  +  vt  (A.l) 

where  all  quantities  refer  to  time  t ,  and  yt  is  the  vector  of  differences  between  the 
observations  and  their  theoretical  values  calculated  from  a  priori  values  of  the  param¬ 
eters;  x 1  is  the  vector  of  adjustments  to  the  a  priori  values  of  the  parameters;  At  is  the 
matrix  of  partial  derivatives  which  relate  the  change  in  parameter  values  to  changes 
in  the  values  of  the  observables;  and  vt  is  a  vector  of  residuals  which  represent  the 
noise  in  the  observations.  The  dynamics  of  the  parameters  are  represented  by  the  state 
transition  equation: 

xt+l  =  Stxt  +  Wt  {A.  2) 

where  xt+i  is  the  vector  of  values  of  the  parameters  at  time  t  + 1;  St  is  the  state 
transition  matrix  at  time  t ,  i.e.,  the  matrix  which  gives  the  expected  changes  in  the 
parameters  between  times  t-f-1  and  t;  and  wt  is  the  vector  of  random  perturbations  to 
the  state.  The  elements  of  wt  corresponding  to  non-stochastic  parameters  are  zero, 
i.e.,  there  are  no  perturbations  of  the  state  with  time.  These  parameters  typically,  but 
not  necessarily,  include  the  positions  of  the  radio  telescopes  and  the  radio  sources,  and 
the  orientation  of  the  terrestrial  coordinate  system  in  space.  The  stochastic  parameters 
typically  include  the  components  of  the  random  walks  and  integrated  random  walks 
being  used  to  represent  the  fluctuations  of  the  clocks  and  the  atmospheric  delays. 

Equations  (A.l)  and  (A. 2)  represent  a  general  set  of  equations  for  linear  dy¬ 
namic  system.  In  order  to  form  the  Kalman  filter  estimator  we  make  the  following 
assumptions: 


<V|)  =  o, 

for  all  t 

(A.  3) 

(*i®i+j)  =  o, 

for  all  j  /0 

(A.4) 

(vtwt+j)  =  0, 

for  all  t  and  j 

(A.5) 

(vtxt+>)  =  0, 

for  all  t  and  j 

(A.6) 

{wt)  =  o, 

for  all  t 

(A.7) 

(wtit’t+i)  =  0, 

for  all  j  -A  0 

(.4.8) 

o' 
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for  all  j  >  0 

(A.9) 
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We  also  define 


(A.l  0) 
(j4.11) 


(vtvt)  =  Vt 

(wtwt)  =  Wt 

The  assumptions  (A. 5)  and  (A.6)  ensures  that  the  measurement  process  and  the  ran¬ 
dom  motion  of  the  system  are  uncorrelated.  The  assumptions  (A. 8)  and  (A.9)  ensure 
that  the  current  state  of  the  system  do  not  affect  the  random  perturbations  in  the 
system  at  later  times.  For  the  classes  of  stochastic  processes  discussed  in  Appendix  B 
this  assumption  is  satisfied. 

The  Kalman  filter  estimation  is  carried  out  sequentially.  The  sequence  for  adding 
the  observations  at  time  1+1  to  the  estimation,  given  the  state,  x‘,  and  its  covariance 
matrix,  Cj,  at  time  1  is, 

Prediction: 

*t+i  —  StxJ 

c;+i  =  s<c'stT^w, 

Update: 

=  *l+i  +  K(yt+i  -  AtTx^+1)  (A. 14) 

Cltl  =  Ct*+1  -  KAt  +  ,C;+1  (A. 15) 

where  K,  the  Kalman  gain,  is  given  by 

K  =  Ci+1AtT+x(Vt+1  +  At+1C»+1Af+1)"1  (A.16) 

Equations  (A.12)-(A.16)  form  a  complete  sequence  for  adding  observations  into  the 
Slter.  When  the  computations  at  time  1+1  are  completed,  the  sequence  is  repeated 
with  quantities  at  time  1  + 2  substituted  for  those  at  time  1+1,  and  those  at  times 
1+1  for  those  at  time  1.  The  sequence  is  repeated  until  all  observations  have  been 
included.  To  start  the  filter,  a  priori  values  for  the  parameters  and  their  covariance 
matrix  are  used  for  quantities  with  sub-  and  superscript  1.  We  refer  to  this  sequence 
as  the  forward  Kalman  filter. 


(A.12) 

(A.13) 
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After  all  observations  have  been  added  to  the  forward  Kalman  filter,  the  final  es¬ 
timates  of  the  non-stochastic  parameters  are  determined.  To  determine  the  estimates 
of  the  stochastic  parameters,  a  backward  or  smoothing  Kalman  filter  solution  must  be 
performed.  This  is  necessary  because  the  estimates  of  stochastic  parameter  for  all  but 
the  last  epoch  of  observations  do  not  contain  the  information  about  the  parameters  es¬ 
timates  supplied  by  later  observations.  Algorithms  for  the  backward  Kalman  filter  do 
not  often  appear  in  the  literature.  (This  is  because  Kalman  filters  are  most  often  used 
in  real  time  applications,  and  only  the  values  of  the  parameters  at  the  current  time  of 
interest.  Liebelt  (1967,  Chapter  6)  gives  an  algorithm  for  the  backward  Kalman  filter 
which  requires  relatively  large  amount  of  computer  memory  to  implement.  We  have 
instead  implemented  the  algorithm  of  Gelb  (1974,  Chapter  5). 

The  Gelb  algorithm  for  smoothing  takes  the  mean  of  the  estimates  of  the  pa¬ 
rameters  from  the  forward  filter  and  from  a  forward  filter  with  time  running  in  reverse 
order.  This  mean  is  taken  at  an  appropriate  position  in  the  filter  equations  such 
that  when  the  mean  is  computed,  all  available  information  about  the  parameters  is 
included.  The  mean  is  taken  using  the  full  covariance  matrices  for  the  forward  and 
backward  estimates  of  the  parameters.  The  choice  of  exactly  which  positions  in  the 
filter  equations  are  used  to  obtain  the  parameter  estimates  for  taking  the  mean  is 
arbitary  provided  that  the  condition  that  all  information  is  used  is  satisfied.  We  have 
chosen  the  positions  corresponding  to  Equations  (A. 14)  and  (A. 15)  for  the  forward 
running  filter,  and  Equations  (A. 12)  and  (A. 13)  for  the  backward  running  filter.  If 
a  back  solution  is  to  be  run  (which  is  not  always  the  case),  then  during  the  forward 
solution  the  parameter  estimates  and  their  covariance  matrix  at  time  t  + 1  are  saved  in 
a  disk  file.  During  the  backward  running  filter,  the  smoothed  estimates  of  the  param¬ 
eters  are  obtained  after  the  prediction  step  from  time  t+2  to  time  t+1.  (Note:  that 
since  the  time  is  running  backwards  in  this  filter  the  observations  from  time  t+2  are 
processed  before  those  from  time  t  +  1.)  The  equations  we  use  for  taking  the  mean  can 
be  derived  from  the  Kalman  filter  equations  given  above.  They  are 

B  =  C+(C^  +  C  +  )"1 
x„  =  x+  +  B(x_  -  x+) 

CB  =  C+  -BC+ 
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(A.17) 

(A.18) 

(A.19) 


where  B  is  the  Kalman  filter  gain  for  taking  the  mean;  x+  and  C*  are  x\  and  Cj  from 
the  forward  running  filter;  x_  and  C_  are  xj+1  and  Ct  +  1  from  the  backward  running 
filter;  xs  and  CB  are  the  smoothed  estimates  of  the  parameters  and  their  covariance 
matrix  at  time  t. 

To  implement  the  smoothing  filter  rigorously  the  Kalman  filter  equations  should 
should  be  written  in  a  more  symmetrical  form  than  that  given  by  Equations  (A. 12)- 
(A.16),  i.e.,  an  a  posterori  contraint  should  be  placed  on  the  parameter  estimates 
after  all  of  the  observations  have  been  included.  The  a  posterori  constraint  would 
then  be  used  as  the  a  priori  constraint  when  the  backward  running  filter  is  started. 
We  have  not  implemented  the  smoothing  equations  in  this  manner.  Instead,  we  place 
very  weak  constraints  on  the  parameters  when  the  backward  running  filter  is  started. 
These  constraints  should  be  weak  enough  not  to  affect  the  parameter  estimates  but 
no  so  weak  as  to  cause  significant  rounding  error  in  the  computations.  The  constraint 
which  we  have  found  works  well  in  practice  is  to  start  the  backward  filter  with  apriori 
variances  which  are  one  thousand  times  greater  than  the  variances  of  the  estimates 
of  the  parameters  at  the  end  of  the  forward  filter  run.  This  weak  constraint  should 
change  the  smoothed  estimates  of  the  parameters  by  about  one  thousandth  of  their 
standard  deviation  (see  discussion  below). 

The  Kalman  filter  equations  given  in  Equations  (A.l)-(A.18)  are  relatively  easy 
to  implement  in  a  computer  program.  However,  there  are  items  about,  and  some 
special  forms  of,  these  equations  which  are  worth  noting.  The  two  items  which  we 
will  discuss  are  the  effects  of  a  priori  constraints  placed  on  the  parameters,  and  the 
implementation  of  white  noise  stochastic  process.  The  Kalman  filter  equations  are 
formulated  such  that  the  covariance  matrix  of  the  parameters  is  manipulated  as  ob¬ 
servations  are  input  to  the  filter.  Consequently,  we  need  to  start  the  filter  with  the 
covariance  matrix  of  the  a  priori  estimates  of  the  parameters.  In  many  applications 
we  would  like  to  start  the  filter  with  no  constraints  placed  on  the  parameters  at  all. 
However,  in  a  computer  this  option  is  not  possible,  i.e.,  infinity  cannot  be  truly  rep¬ 
resented.  In  addition,  if  we  look  at  Equation  (A. 15),  placing  large  a  priori  constraints 
will  introduce  large  rounding  errors  into  the  processing  since  the  final  covariance  ma¬ 
trix  is  obtained  by  decrementing  as  each  new  set  of  observations  is  added  to  the 

filter.  Given  these  considerations,  it  is  worthwhile  to  consider  how  the  a  priori  covari¬ 
ances  affect  the  final  estimates  of  the  parameters  and  their  covariance  matrix.  Since 
any  a  priori  constraint  can  be  treated  as  a  set  of  observations,  we  will  use  the  Kalman 


filter  equations  to  investigate  the  effects  of  the  constraints.  For  the  constraints,  the 
partial  derivative  matrix  At  is  simply  a  unit  matrix,  and  the  Kalman  gain  when  the 
constraints  are  applied  will  be  given  by, 

K  =  C^Co  +  C*)-1  ( A.20 ) 

where  Co  is  the  matrix  of  a  priori  variances,  and  C*  is  the  covariance  matrix  of  the 
parameter  estimates  without  the  constraint  applied.  If  we  use  the  matrix  identity  that 
(A  +  B)-1  =  B-1(AB_1  +  1)  l,  where  I  is  a  unit  matrix,  then  K  reduces  to, 

K  =  (C0(C*)_1  +  I)1  (A. 21) 

and  the  estimates  of  the  parameters  with  the  constraint  applied  becomes 

x  =  x\  +  (Co(Ct)-1  +  l}~\x0  -  x\)  (A. 22) 

where  Xo  is  a  vector  of  u  priori  values.  In  general  Equation  (A. 22)  is  difficult  to  solve 
analytically.  However,  we  can  treat  the  special  case  when  C*  is  a  diagonal  matrix. 
(We  also  assume  that  C0  is  a  diagonal  matrix.)  Denoting  the  diagonal  elements  of  Co 
and  Cj  for  the  tth  parameter  by  cTq,  and  o^,  respectively,  the  change  in  the  estimate 
of  the  ith  parmeter  when  the  constraint  is  applied  is  given  by 

Ax,  =  — ^-(xoj  —  Xft)  (A. 23) 

°o 

pi 

where  we  have  assumed  that  i  «  1.  From  Equation  (A. 23)  we  see  that  the  change 
in  the  parameter  estimate  due  to  the  constraint  being  applied  is  the  adjustment  to 
parameter  value  by  the  ratio  of  the  a  posterori  variance  to  the  a  priori  variance. 
Similarly,  it  may  be  shown  the  change  in  the  estimate  of  the  variance  of  the  parameters 
is  also  given  by  -A.  Another  case  which  can  be  easily  studied  is  when  all  of  the 

°Q 

parameters  are  highly  correlated.  In  this  case,  the  effect  of  the  constraint  is  n  times 
greater  than  in  the  uncorrelated  case,  where  n  is  the  number  of  highly  correlated 
parameters.  We  have  found  that  these  simple  rules  work  well  for  any  Cj  encountered 
in  our  processing. 

The  special  form  of  the  Kalman  filter  equations  which  we  wish  to  study  is  the 
case  of  white  noise  stochastic  processes.  For  this  type  of  process  the  part  of  the  state 
transition  matrix  for  the  white  noise  process  is  zero.  We  will  also  assume  the  white 
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noise  process  is  uncorrelated  with  any  of  the  other  noise  in  the  system.  We  will  now 
rewrite  the  Kalman  filter  equations  in  partitioned  form.  Equation  (A.l)  becomes 


yt  -  (At  A't )  ^  +  vt 


(A.  24) 


where  n*  is  the  vector  of  parameters  in  the  white  noise  process,  and  A\  is  the  matrix 
relating  nt  to  yt.  The  other  quantities  in  Equation  (A. 23)  have  the  same  meaning  as 
before.  The  state  transition  Equation  (A. 2)  becomes 


+  _  /St  °Wit\  (  wt  \ 

\nt+ \)  V  0  0/  \Mt/ 


(A. 25) 


The  application  of  the  Kalman  filter  equations  to  this  type  of  system  is  straight  for¬ 
ward.  The  prediction  step  for  the  parameters  and  their  covariance  matrix  is  given  by 
(from  Equations  (A. 12)  and  (A. 13)) 


r,t 

v-n+i 


Xt+  1 

”t  +  l 

0 

1 

^  t+1 


‘t 


n| 

0 

0 


■( 


St  0 
0  0 

StC{S tT 
0 

stc\s?  +  wt 

0 


(A.26) 


+ 


Wt 

0 


0 

Nt 


0 

Nt 


(A.  27) 


where  C't+1  is  the  submatrix  of  covariance  matrix  for  the  white  noise  process.  It 
should  be  noticed  that  this  submatrix  is  simply  Nt.  The  Kalman  gain  for  this  system 
reduces  to 

(*,)  =  (  )  (Q.  +  a'.NiA',T  +  A,c‘+1Af )"  (/t.28) 


where  K'  is  the  submatrix  of  the  Kalman  filter  gain  for  the  white  noise  process.  If  we 
°xamine  the  expression  for  the  Kalman  gain  of  the  non-white  noise  process  parameters 
in  Equation  (A.28)  we  notice  that  the  same  expression  for  the  Kalman  gain  would  have 
been  generated  had  we  not  explicitly  included  the  white  noise  process  parameters  in 
the  filter,  but  had  instead  simply  modified  the  covariance  matrix  of  the  observations  to 
include  the  effects  of  the  white  noise  process  «.e.  Qt  + A  tNtA't  would  be  substituted 
for  Qt.  This  modification  is  valid  for  the  remaining  equations  in  the  Kalman  filter  as 
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well.  We  have  found  this  procedure  useful  whenever  white  noise  processes  are  needed  in 
the  Kalman  filter  since  it  reduces  the  number  of  parameters  which  need  to  be  explicitly 
included  in  the  filter. 

A  very  useful  application  for  the  above  procedure  occurs  when  both  delay  and 
rate  data  are  included  in  the  filter,  and  a  random  walk  (for  delays)  is  included  in  the 
stochastic  models.  For  the  rate  data,  this  random  walk  process  will  manifest  itself 
as  white  noise.  If  the  duration  of  the  observations  is  short  compared  to  time  interval 
between  the  observations,  then  the  white  noise  at  the  “instant”  of  observation  will 
be  independent  of  the  integrated  effect  of  the  white  noise  sensed  by  the  delay  data. 
Thus,  the  effects  of  the  random  walk  on  the  rate  data  can  be  included  using  the 
procedure  given  above.  However,  in  practice,  the  duration  of  the  observations  is  not 
short  compared  to  interval  time  between  measurements  and  thus  the  procedure  given 
about  is  not  valid.  We  will  defer  furthur  discussion  of  this  topic  until  we  discuss  the 
effects  of  the  finite  of  duration  of  observations  on  modeling  stochastic  processes  in 
Appendix  B. 


Appendix  B.  Properties  of  stochastic  processes 

Here  we  discuss  the  properties  of  the  stochastic  processes  we  have  used  in  out 
implementation  of  the  Kalman  filter.  We  will  firstly  define  each  of  the  statistical 
characterizations  given  in  Table  A.l;  then  we  will  discuss  the  derivations  of  some  of 
these  characterizations;  and  finally  we  will  discuss  the  effects  of  averaging  the  stochastic 
processes.  For  each  of  the  stochastic  processes  we  give  the  following,  (i)  Driving 
differential  equation.  In  all  cases  the  process  is  assumed  to  be  driven  by  white  noise, 
w{t).  (ii)  Solution.  This  expression  is  the  solution  to  the  differential  equation  which 
drives  the  process.  (Hi)  covariance  function,  R(t\,t 2).  The  covariance  function  is 
defined  as, 

R{tuh)  =  (p[ti)p{t2))  [B.  1) 

where  p(t)  is  the  value  of  the  process  at  time  t,  and  ()  denotes  the  expectation  value. 
When  the  covariance  function  depends  only  on  ti~t2,  the  function  is  stationary,  and 
the  covariance  function  is  expressed  in  terms  of  r=tj-t 2.  (iv)  Power  spectral  density 
function  (PSD),  S(f).  For  stationary  stochastic  processes,  the  PSD  is  defined  as, 

S(f)  =  /  R(r)e-,2KfTdr  {B .  2) 

J  OO 
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where  /  is  frequency  in  cycles  per  unit  time,  (v)  2-sample  Allan  variance,  o^(t).  The 
2-sample  Allan  variance  is  defined  by  (Vessot,  1976), 

tf(T)  =  <P2(<  +  T)  ~  2p(<  +  r)p(t)  +  p2{t))/r 2  {B. 3) 

where  the  expectation  is  taken  over  all  times  t.  For  all  of  the  stochastic  processes  given 
in  Table  A.l,  <7y(r)  is  only  a  function  of  r,  the  time  interval  between  samples,  and  not 
the  time  at  which  the  samples  are  taken.  The  form  of  the  2-sample  Allan  variance  is 
such  that  the  results  are  independent  of  any  mean  value  or  mean  slope  in  the  values 
of  p{t).  This  property  is  particularly  advantageous  to  the  study  of  hydrogen  maser 
clocks  because  their  absolute  frequency  is  often  not  well  controlled.  Also,  we  refer  to 
the  Allan  standard  deviation  most  of  the  time  which  we  define  to  be  the  square-root  of 
the  2-sample  Allan  variance,  (vi)  Structure  function,  D(t,r).  The  structure  function, 
which  is  often  used  in  atmospheric  studies,  is  defined  as 

C((,r)  =  ((p(I  +  r)-p(<))1).  (fl.4) 

where  the  expectation  is  taken  over  all  times  t.  As  for  the  Allan  variance,  the  structure 
function  may  depend  on  t  as  well  as  r. 

The  derivations  of  the  entries  given  in  Table  A.l  are  tedious  but  straightforward 
given  the  solutions  to  the  differential  equations  for  each  stochastic  process.  For  each 
of  the  processes,  $Q  (where  a  is  w  for  the  white  noise,  r  for  the  random  walk,  t  for 
the  integrated  random  walk,  and  g  for  the  fint-order  Gauss  Markov  processes)  is  the 
PSD  of  the  white  noise  driving  the  process.  The  variance  of  the  driving  white  noise  is 
$0,6(0)  where  6(0)  is  the  Dirac  delta  function  i.e.,  (w(t)w(t  +  r))  =  $a6(r). 

Some  general  comments  can  be  made  about  the  nature  of  the  processes  given  in 
Table  A.l.  For  the  random  walk  and  the  integrated  random  walk  the  PSD  functions 
are  not  defined  because  neither  of  these  processes  is  stationary,  i.e.,  their  variances 
tend  to  infinity  as  time  tends  to  infinity.  In  our  application  this  property  does  not  pose 
a  major  problem  because  we  only  need  the  stochastic  model  to  predict  the  variance  of 
the  change  in  the  process  over  the  interval  of  time  between  observations.  An  interval 
which  is  typically  less  than  1000  seconds.  When  there  are  large  gaps  in  the  VLBI 
data  (>1000  seconds),  the  random  walk  and  integrated  random  walk  models  tend 
to  overestimate  the  variance  of  the  changes  in  the  processes.  The  net  effect  is  that 


information  about  the  change  in  the  process  values  over  the  gap  will  come  more  from 
the  VLBI  data  themselves  than  from  the  prediction  based  on  the  stochastic  model. 

In  our  implementation  of  the  Kalman  filter  the  random  walk  process  is  the  one 
most  often  used.  This  process  seems  appropriate  for  representing  the  short  period 
fluctuations  (<  104  seconds)  of  the  clocks  and  the  fluctuations  of  the  atmospheric 
delays.  Since  the  changes  in  the  values  of  the  random  walk  processes  over  the  duration 
of  a  VLBI  observation  are  not  small,  we  will  consider  in  detail  the  incorporation  of 
this  stochastic  process  in  the  filter.  We  consider  two  VLBI  observations  of  durations 
Ai  and  A2  with  the  centers  of  the  observations  at  times  and  £2-  The  estimates  of 
the  values  of  VLBI  delay  and  rate  observables  averaged  over  the  two  durations  is 


,  =  /  p(u)du/2A, 

Jt,- A. 

3  ft,+A, 

*  =  2A3  Jt  A  ^  ~  *‘)p(u)du 


(B.5) 
(B.  6) 


where  t  is  either  1  or  2  denoting  the  first  or  second  observation.  Given  these  expressions 
for  the  delay  and  rate,  we  can  subsitute  the  equations  for  the  random  walk  process 
noise  (from  Table  A.l),  and  evaluate  the  statistical  properties  of  the  averages.  The 
following  expressions  can  be  derived  for  the  statistical  properties  of  the  delays  and 
rates  at  the  two  times: 


(r,2)  =  4>rt,  ~  $rA,/6 

[B.  7) 

<r2)  =  6*r 

5A, 

[B.  8) 

{UU)  =  $r/2 

(B- 9) 

(t»T2>  =  for  *i  <  t2) 

(5.10) 

(fir2)  =  0 

(5.11) 

(tlT2)  =  $r. 

(5.12) 

Our  implementation  of  the  Kalman  filter  docs  not  adhere  exactly  with  these  equa¬ 
tions.  We  neglect  the  small  decrease  in  the  variance  of  (r2)  due  to  the  term  $rA,/6 
because  the  appropriate  durations  of  the  observations  is  not  always  known  especially 


for  the  VLBI  data  taken  in  the  1970’s.  Nor  do  we  account  for  the  small  correla¬ 
tion  (5A,/y/l2(t,  -  A,/6),  t,  >  A,/2)  between  the  delay  and  rate  observations  at 
each  epoch  of  observation.  To  rigorously  include  the  effects  of  this  correlation  in  the 
Kalman  filter  would  require  the  introduction  of  an  additional  parameter  for  each  ran¬ 
dom  walk  process.  This  parameter  would  represent  the  average  rate  of  change  of  the 
random  walk  over  the  duration  of  an  observation.  The  process  would  still  be  a  white 
noise  sequence,  but  the  sequence  would  be  correlated  with  the  random  walk  sequence 
for  the  delays.  Thus,  the  assumption  used  in  deriving  Equation  (A. 27)  would  not 
valid,  and  we  could  not  simply  incorporate  the  rate  process  into  the  rate  data  co- 
variance  matrix.  We  do  not  feel  that  this  is  a  major  problem  in  the  Kalman  filter 
implementation  because  we  are  neglecting  a  small  amount  of  information  about  the 
rate  data  which  should  result  in  the  Kalman  filter  yielding  slightly  larger  variances 
for  the  estimated  parameters.  If  in  the  future,  the  time  interval  between  observations 
becomes  very  short  compared  to  the  duration  of  the  observations,  we  may  modify  the 
Kalman  filter  implementation  to  incorporate  the  correlation  between  delay  and  rates. 
Also,  at  this  time,  we  generate  all  of  our  geodetic  results  using  only  delay  data.  The 
neglected  correlation  has  no  impact  at  all  on  this  type  of  solution.  We  adopt  this 
procedure  because  including  the  rate  into  the  geodetic  solutions  has  a  minimal  effect 
on  the  estimates  of  the  geodetic  parameters  and  their  variances. 


. 

TABLE  A.l.  Characteristics  of  some  stochastic  processes. 


White  noise,  pw(t) 

Random  Walk,  pr(t) 

Pw{t)=  u;(t) 

pw[t  4-  At)—  w(t  4-  At) 

R{t)=  $w6(r) 

S(f)=  *w 

cr£(r)=  3<t>w/r2 

D{r)  =  2$w6(0) 

pw(t  -f  At)=  w(t  4-  At) 

At 

pr(t  4-  At)=  pr(f)  4-  /  w(u  +  t)du 

0 

R{ti,t2)=  mm(t1,t2) 

S(f)  not  defined 

*y(T)=  *r/r 

D(t)=  $rT 

At 

pr(t  4- At)=  pr(t)  4-  /  w(u  +  t)du 

0 

Integrated  Random  Walk,  p,(t) 

First  order  Gauss  Markov,  pg(t) 

Pi{t  4-  At)—  p,(t)  4-  Atp,(t) 

At 

4-  f  (At  -  u)w(u  4-  t)du 

0 

Pi(t  4-  At)  =  p,(t )  4-  /0A<  to(u  4-  t)du 
R{tut2)=  (3t2-t1)/6 

S(f)  not  defined 

°1{T)=  *iT/ 3 

D(t,r)=^(r3/3  +  T2t) 

Pi{t  4-  At)=  p,(t)  4-  Atpi(t) 

At 

4-  /  (At  -  u)w{u  4-  t)du 

0 

p,(t  +  A t)=  pi(t)  4-  J0Af  w(u  4-  t)du 

Pg[t  +  At)=  e~At0Pg{t) 

4 -e~At(*  f  eu^w(t  +  u)du 

0 

*w-  It'-™ 

s(n ~  W(i  +(2*f/0i‘) 

=  -  4e-W  + '-W>) 

D{t)~  ^(1-e-l^) 
pff(f  +  At)  =  e~Atppg{t) 

At 

+e~At^  j  e«/3u,((  4.  u)du 

0 

In  each  of  the  boxes  we  give  the  characteristics  of  the  stochastic  processes  discussed  in  the  text. 
The  first  item  in  each  box  is  the  driving  differential  equation  for  the  process.  The  second  entry 
is  the  solution  to  the  differential  equation  for  an  arbitrary  excitation  function.  i?(r)  or  (t  i ,  1 2 ) 
is  the  covariance  function  of  the  process,  assuming  that  the  excitation  of  function,  w{t),  is  white 
noise  with  a  power-spectral  density  (PSD)  of  where  q  is  w  for  the  white  noise  process,  r 
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for  the  random  walk,  i  for  the  integrated  random  walk,  and  g  for  the  first-order  Gauss  Markov 
process.  S[f)  is  the  PSD  function  where  /  is  in  cycles  per  unit  time.  ( r )  is  the  2-sample 
Allan  variance  at  sampling  time  r.  D(t,r)  or  D(t)  is  the  structure  function  for  time  interval  r. 
The  final  entry  for  each  process  is  the  solution  to  the  driving  differential  equation  for  arbitrary- 
excitation  function  w(f).  The  definitions  for  each  of  these  quantities  and  some  discussion  of  the 
entries  given  in  the  table  are  given  in  Appendix  B. 
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METHODS  OF  CORRECTION  FOR  THE  “WET"  ATMOSPHERE  IN  ESTIMAT¬ 
ING  BASELINE  LENGTHS  FROM  VLBI 


G.  Elgered,  J.  L.  Davis,  T.  A.  Herring,  and  I.  I.  Shapiro 
Center  for  Astrophysics 
60  Garden  Street 
Cambridge,  MA  02138 


ABSTRACT.  The  error  in  VLBI  estimates  of  baseline  length  caused  by  unmodelled  variations 
in  the  propagation  path  through  the  atmosphere  is  greater  for  longer  baselines.  We  present  and 
discuss  series  of  estimates  of  baseline  lengths  obtained  using  different  methods  to  correct  for  the 
propagation  delay  caused  by  atmospheric  water  vapor.  The  main  methods  are  use  of  data  from 
a  water-vapor  radiometer  (WVR)  and  Kalman-filtering  of  the  VLBI  data  themselves  to  estimate 
the  propagation  delay.  Since  the  longest  timespan  of  WVR  data  associated  with  geodetic  VLBI 
experiments  was  obtained  at  the  Onsala  Space  Observatory  in  Sweden,  we  present  results  for  the 
following  three  baselines:  (l)  Onsala-Wettiell,  FRG  (920  km),  (2)  Onsala-Haystack/Westford, 
MA  (5600  km),  and  (3)  Onsala-Owens  Valley  (7914  km). 

1.  INTRODUCTION 

This  paper  addresses  the  utility  of  water-vapor  radiometers  (WVR)  in  geodetic  radio- interferometry 
experiments.  The  WVR  data  are  providing  a  prton  information  about  the  wet  delay.  It  is  also 
possible  to  estimate  a  correction  to  any  a  prion  delay  during  post-processing.  The  estimate  can  be 
of  a  mean  senith  bias  for  the  entire  experiment  or  of  values  of  samples  of,  say  an  assumed  random 
(Markov)  process.  We  have  analysed  77  experiments  (made  during  1980-1985)  several  times,  each 
time  with  a  different  method  to  correct  for  the  wet  delay  at  Onsala,  but  with  the  atmospheric 
delays  for  all  other  sites  (as  well  as  the  clocks)  modeled  as  Markov  processes. 


2.  RESULTS  AND  DISCUSSION 

The  results  are  presented  in  Table  I  as  weighted  root-mean -square  (WRMS)  scatters  of  baseline 
lengths  about  estimated  slopes.  The  WRMS  value  is  a  measure  of  repeatability.  We  also  present 
one  solution  where  observations  made  at  elevation  angles  lower  then  25°  (c  <  25°)  at  Onsala  were 
deleted  since  low  elevation  angle  observations  are  not  important  when  no  delay  corrections  are 
estimated  ( Htrring  1986).  The  WRMS  for  the  Wettsell— Onsala  baseline  is  given  with  respect  to 
its  mean  value  since  no  baseline  change  is  expected.  For  this  baseline,  when  the  WVR  data  are 
used  in  place  of  the  Markov  process,  the  WRMS  decreases  from  5.2  to  4.0  mm.  If  we  assume  that 
atmospheric  delay  errors  are  not  con-elated  with  other  enors,  and  that  the  wet  atmospheres  over 
Wettsell  and  Onsala  are  statistically  similar,  the  WRMS  would  be  2.2  mm  were  a  WVR  installed 
at  Wettsell.  The  results  for  the  Haystack/Westford-Onsala  baseline  are  presented  for  each  of 
the  former  antennas  separately  since  the  WRMS  is  significantly  different  in  the  two  cases.  This 
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TABLE  I.  BASELINE  LENGTH  REPEATABILITY 


Baseline 

Onsala 

to 

Method  used  to  correct  for 
the  wet  delay  at  Onsala 

A  priori  Adjustment 

Mean  baseline3 
+919660  m 
(mm) 

WRMi 
about  mean 
(mm) 

WettzeU1 

None 

Offset 

999  ±  1 

5.8 

WettzeU1 

None 

Markov 

1000  ±  1 

5.2 

WettzeU1 

WVR 

Offset 

999  m  1 

5.1 

WettzeU1 

WVR 

None 

996  ±  1 

4.0 

WettzeU1 

WVR  (c  >  25s) 

None 

996  ±  1 

4  4 

Slope5 

WRMS  about  slope 

(mm  /year) 

(mm) 

Haystack3 

None 

Offset 

17.8  ±  1.4 

14  1 

Haystack3 

None 

Markov 

19.0  ±  1.2 

12.8 

Haystack3 

WVR 

Offset 

16.1  ±  1.2 

12.5 

Haystack3 

WVR 

None 

14.2  ±  1.5 

15.0 

Haystack3 

WVR  (c  >  25") 

Non* 

16.6  ±  1.1 

11.5 

Westford3 

Non* 

Offset 

15.9-2.5 

21.6 

Westford3 

None 

Markov 

17.8m  2.1 

18.3 

Westford3 

WVR 

Offset 

14.5m  2.5 

21.8 

Westford3 

WVR 

None 

13.7-2.8 

23.1 

Westford3 

WVR  (c  >  25  =  ) 

None 

15.8  ±2.3 

18.6 

Owens  VaUey4 

None 

Offset 

12.8  ±4.0 

39.1 

Owens  VaUey4 

Non* 

Markov 

12.7  ±  3.2 

31  4 

Owens  Valley4 

WVR 

Offset 

10.3  ±  3.2 

31.3 

Owens  VaUey4 

WVR 

None 

90:6.1 

57.6 

Owens  VaUey4 

WVR  (e  >  25") 

None 

11.4  ±  3.6 

35  4 

1  25  experiments  3  36  experiments  3  45  experiment*  4  29  experiment* 

6  The  sigmas  we  scaled  eo  that  reduced  x2  *=  1  ( Herring  ct  at.  1986). 


difference  indicates  that  the  accuracy  of  the  Westford-Onsala  baseline  estimates  is  not  Limited  by 
atmospheric  errors.  It  is  also  dear  that  low  elevation  observations  do  not  improve  the  accuracy 
of  baseline  length  estimates  when  WVR  data  are  used  and  no  adjustment  to  the  a  pnon  delay  is 
estimated.  Finally,  for  the  longest  baseline  (Owens  Valley-Onsala)  it  appears  better  to  estimate 
a  correction  to  the  aenith  delay  inferred  from  the  WVR  data  rather  than  to  discard  low  elevation 
observations,  because  a  bias  affects  the  accuracy  of  the  vertical  position  which  is  more  important 
for  estimates  of  length  of  longer  baselines. 


3.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  WORK 


The  differences  in  the  rates  obtained  with  the  different  methods  are  small  compared  to  1  cm/year. 
It  is  important  to  minimize  possible  biases  in  the  atmospheric  delays  inferred  from  WVR  data 
A  main  source  of  such  bias  is  the  uncertainty  of  the  attenuation  coefficients  due  to  water  vapor. 
Another  important  task  is  to  find  the  'best*  elevation  cut-off  angle  when  WVR  data  are  to  be 
used. 
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ABSTRACT 


An  important  source  of  error  in  very-long-baseline  interferometry  (VLB1) 
estimates  of  baseline  length  is  unmodeled  variations  of  the  refractivity  along  the 
propagation  path  through  the  atmosphere.  We  present  and  discuss  the  method 
of  using  data  from  a  water-vapor  radiometer  (WVR)  to  correct  for  the  propaga¬ 
tion  delay  caused  by  atmospheric  water  vapor.  Data  from  different  WYR's  are 
compared  with  estimated  propagation  delays  obtained  by  Kalman-filtering  of  the 
VLBI  data  themselves.  The  consequenses  of  using  either  WVR  data  or  Kalman 
filtering  to  correct  for  atmospheric  delays  at  the  Onsala  VLBI  site  are  investigated 
by  studying  the  repeatability  of  estimated  baseline  lengths  from  Onsala  to  several 
other  sites.  The  lengths  of  the  baselines  ranges  from  91S  to  7941  km.  The  repeata¬ 
bility  obtained  for  baseline-length  estimates  show  that  the  methods  of  water- vapor 
radiometry  and  Kalman  filtering  offer  comparable  accuracies  when  ap-plied  to  the 
climate  of  the  Swedish  west  coast.  It  is  also  clear  that  observations  made  at  low 
elevation  angles  should  not  be  used  at  a  site  where  WVR  data  are  available  if  no 
estimations  of  the  atmospheric  delay  is  made  for  that  site  from  the  interferometer 
data.  The  actual  “best”  cut-off  in  elevation  angle  depends  on  the  accuracy  of  the 
total  delay  since  its  error  increase  with  increasing  air  mass.  In  this  study  the  best 
cut-off  is  found  to  be  about  20°  when  WVR  data  are  used. 


INTRODUCTION 


To  correct  for  the  “wet”  atmospheric  delay,  water- vapor  radiometers  (WYR's) 
have  been  developed  for  radio-interferometry  experiments  designed  for  estimation 
of  geodetic  parameters.  This  delay,  although  much  smaller  then  the  delay  caused 
by  the  “dry”  constituents  of  air,  is — due  to  its  variability  in  space  and  time — a 
major  source  of  error  in  estimates  of  geodetic  parameters  such  as  baseline-lengths. 

WVR  data  can  be  treated  as  a  priori  ( i.e .,  prior  to  least-squares  estimation 
of  site  positions,  etc.)  information  about  the  wet  delay.  Other  a  priori  data  are 
ground  measurements  of  humidity  and  temperature  that  can  be  used  with  a  model 
to  predict  the  wet  delay.  Due  to  poor  mixing  of  wet  and  dry  air  the  accuracy  of 
this  type  of  model  is  expected  to  be  too  poor  to  be  useful  for  our  application.  For 
comparison,  however,  we  use  a  model  of  this  type:  Saasiamoinen't  [19721;  we  will 
refer  to  it  as  the  ground-based  model. 

It  appears  that  there  is  no  practical  possibility  other  than  use  of  a  remote 
sensing  instrument,  such  as  the  WVR.  in  connection  with  baseline-length  esti¬ 
mation  from  very  long-baseline-interferometry  (VLBI)  data.  However,  it  is  also 
possible  to  estimate  a  correction  to  any  a  priori  delay  together  with  geodetic  and 
other  parameters  from  the  VLBI  data.  We  can  estimate  a  mean  zenith  bias  for 
the  entire  experiment  or  values  of  samples  of,  say,  an  assumed  random  (Markov) 
process  [Herring  et  al ,  1988]  representing  the  delay. 

This  paper  addresses  the  utility  of  water- vapor  radiometers  (WVR)  in  geode¬ 
tic  radio-interferometry  experiments.  First,  a  general  background  of  atmospheric 
delays  and  water- vapor  radiometry  will  be  given.  Thereafter,  we  describe  the 
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different  WVR’s  that  have  been  used  to  collect  the  data  now  analyzed.  The  re¬ 
sults  are  presented  in  two  ways.  First  we  compare  the  inferred  wet  delay  variations 
from  the  WYR  and  the  Kalman  filter  for  several  experiments  Second,  we  compare 
the  repeatability  of  estimated  baseline  lengths  obtained  using  different  methods 
to  correct  for  the  atmospheric  delays.  Since  the  longest  timespan  of  WVR  data 
associated  with  geodetic  VLBI  experiments  was  obtained  at  the  Onsala  Space 
Observatory  in  Sweden,  we  present  series  of  estimates  of  lengths  from  Onsala  to 
Wettzell,  FRG;  Haystack/Westford,  MA;  Richmond,  FL;  Owens  Valley,  CA;  and 
Fort  Davis,  TX. 

THE  ATMOSPHERIC  DELAY  AND  WATER  VAPOR  RADI OMETRY 
The  delay  through  the  neutral  atmosphere  depends  on  two  terms  (Davit  et 
o/.,  1985).  The  first  is  called  the  “hydrostatic’'  (or  dry)  delay.  Its  value  AL\  in 
the  zenith  direction,  expressed  in  meter  is 

A Lk  =  (0.0022768  i  0.0000015)  --ft  -  (1) 

/(***)■«) 

where  P0  is  the  total  pressure  at  the  ground  in  mbar  and  where 

/($>,#)  =  (1  -  0.00266 cos 2*  -  0.00028H)  (2) 

is  used  to  model  the  variation  of  the  acceleration  due  to  gravity  with  latitude  $  and 
the  height  of  the  station  H,  in  km,  above  the  geoid.  The  uncertainty  given  for  the 
hydrostatic  delay  is  the  root  sum  square  of  the  uncertainties  in  the  measurements 


of  the  refractivity,  the  acceleration  due  to  gravity,  the  universal  gas  constant, 
and  the  variability  of  the  dry  mean  molar  mass,  but  not  in  any  departures  from 
hydrostatic  equlibrium. 

The  elevation  dependence  of  the  hydrostatic  delay  is  modeled  by  using  a 
“mapping  function”  whose  accuracy  is  improved  if  other  meteorological  measure¬ 
ments,  such  as  of  the  ground  temperatures,  are  used  (see  c.g.,  HopficlJ  [1971]). 
The  mapping  function  used  for  the  hydrostatic  delay  in  this  study  is  presented  by 
Davit  et  al  [1985]. 

The  second  term  of  the  total  delay  is  the  wet  delay,  AL„,  which,  defined 
consistently  with  the  hydrostatic  delay,  can  be  written  as 

[A  oo  ,oc 

17  J  ~ds  +  3.776  x  10Jy  j^ds  (3) 

expressed  in  the  same  units  as  is  the  path,  t ;  T  is  the  temperature  in  K;  and  e  is 
the  partial  pressure  of  water  vapor  in  mbar.  This  wet  delay  is  determined  mainly 
by  the  amount  of  water  vapor  along  the  atmospheric  path. 

The  WVR  measures  the  emission  from  the  sky  at  two  (or  more)  well  sep¬ 
arated  frequency  bands  near  the  water-vapor  emission  line  centered  between  22 
and  23  GHz.  The  WVR  intensity  output,  for  the  band  closest  to  the  center  of  this 
line  will  depend  mainly  on  the  amount  of  water  vapor  in  the  direction  the  WVR 
antenna  is  pointed.  A  second  frequency  band  is  needed  to  correct  for  emission 
caused  by  liquid  water,  which  occasionally  can  be  larger  than  the  vapor  contri¬ 
bution  (even  at  the  center  of  the  water-vapor  line).  Using  the  Rayleigh-Jeans 
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approximation  and  the  equation  of  radiative  transfer  [Chandrasehar.  1S5C\  we 
can  write  the  sky-brightness  temperature  measured  by  the  WYR  as 


T,  =  Ti3e~T~  +  /  T(s)  a(t)e-r(,)ds 

Jo 


(4) 


where  Tj,  is  the  “background”  radiation  (i.e.,  from  outside  the  earth's  atmo¬ 
sphere).  In  the  frequency  band  relevant  to  WVR’s  is  due  to  the  cosmic  back¬ 
ground  radiation.  T  (s)  is  the  physical  temperature  along  the  path  and  a  (*)  is 
the  attenuation  coefficient  due  to  water  vapor,  liquid  water  and  oxygen.  The 
parameter  r  (»)  is  the  corresponding  opacity  from  the  ground  to  the  point  s 


When  $  is  equal  to  infinity  the  opacity  is  written  as  r <*,.  Each  of  the  attenuation 
coefficients  has  its  own  frequency  dependence  which  makes  it  possible  to  separate 
approximately  the  effects  due  to  oxygen,  water  vapor,  and  liquid  water. 

When  studying  (3)-(5)  it  becomes  clear  that  by  making  some  approximations 
we  can  formulate  several  (non-unique)  algorithms  which  will  allow  us  to  estimate 
the  wet  delay  from  WVR  measurements  of  sky-brightness  temperature.  It  is  pos¬ 
sible  to  take  many  different  approaches.  Empirical  parameters  in  the  algorithm 
show  dependence  on  meteorological  conditions  such  as  vertical  profiles  of  temper¬ 
ature  in  the  atmosphere.  Algorithms  can  therefore  be  optimized  to  a  particular 
site.  Different  wet  delay  algo  rithms  have  been  derived  by  Retch  [1983],  Gory  ef 
al  [1985],  Johansson  et  al  [1987],  and  Rohinson  [1988].  These  algorithms  are 
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derived  for  many  \LBI  sites  used  (or  to  be  used)  in  geodetic  experiments  ar.d 
most  of  them  make  use  of  meteorological  parameters  measured  at  the  ground. 

Let  us  define  the  “algorithm  error"  as  the  error  of  the  wet  delay  inferred  from 
noise  free  radiometer  observables.  The  algorithm  error  can  then  be  divided  into 
two  parts.  The  first  part  is  a  scatter  which  varies  on  time  scales  of  hours  or  longer. 
We  obtain  this  error  if  we  assume  that  the  expressions  for  radio  signal  attenuation 
used  to  derive  the  algorithm  are  correct.  For  the  algorithms  referred  to  above, 
the  root-mean-square  (RMS)  error  of  this  part,  for  the  zenith  direction,  varies 
between  0.1  to  0.4  cm,  depending  on  algorithm  and  site,  for  weather  situations 
excluding  rain  or  heavy  rain-clouds.  The  second  part  of  the  algorithm  error,  due 
to  uncertainties  in  the  attenuation  coefficients,  is  believed  to  be  more  of  a  bias 
error  over  time  scales  of  a  day,  the  typical  length  of  a  geodetic  YLB1  experiment. 

The  uncertainty  in  the  attenuation  of  oxygen  is  the  least  important  and 
corresponds  to  an  uncertainty  in  the  inferred  delay  in  the  zenith  direction  of  less 
than  2  mm.  This  value  is  derived  by  comparing  the  results  from  using  different 
attenuation  coefficients  for  oxygen  [Meeks  and  Xi/'ey,  1963;  Snider  and  Westwater. 
1969;  lielc,  1985;  Liebe,  1987;  Rosenkranz.  1987'  in  the  algorithm. 

The  attenuation  due  to  liquid  water  is  assumed  to  be  proportional  to  the 
square  of  the  frequency.  This  frequency  dependence  is  valid  as  long  as  the  size 
of  the  liquid  water  drops  are  much  smaller  than  the  wavelength  of  the  attenuated 
signal  [SlacUn  1966].  For  all  the  WVR’s  used  for  geodetic  YLBI,  the  “cloud  correc¬ 
tion  band”  has  the  highest  frequency  and  is  centered  at  about  31  GHz  (wavelength 
fa  1  cm).  If  the  size  of  the  drops  in  the  atmosphere  become  comparable  to  the 
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wavelength  (say  a  few  tenths  of  a  mm)  the  algorithm  will  overestimate  the  wet 


delay  [  Wcstwaier  1972].  An  example  of  the  effect  cf  large  crop;  is  given  in  Figure 
1,  where  the  effect  is  further  increased  by  the  accumulation  cf  water  drops  on  the 
teflon-covers  of  the  horn  antennas.  During  the  rain,  the  WVR  data  are  then  of  no 
use  for  the  calculation  of  accurate  delay  corrections. 

The  uncertainty  in  the  attenuation  due  to  water-vapor  is  the  most  important 
factor.  The  present  uncertainties  of  these  formulas  are  about  4  to  89c.  with  the 
uncertainty  increasing  with  decreasing  temperature  [ Elgered  1989]  These  uncer¬ 
tainties  are  obtained  by  using  different  expressions  for  the  attenuation  coefficients 
in  the  derivation  of  the  algorithm  [Staehn  19CC.  Waters  197G.  Liele  1985.  Liele 
19871.  Even  a  47c  error  is.  however,  not  negligible.  Zenith  wet  delays  of  between 
10  cm  and  30  cm  are  not  unusual  in  the  temperate  summer.  For  a  typical  elevation 
angles  of  30°,  we  will  then  obtain  errors  of  the  order  of  1-3  cm  in  our  estimated 
delay  corrections.  These  errors  are  large  when  compared  to  the  uncertainties  of 
VLBI  group  delay  observations  which  are  1  cm  or  less. 


BRIEF  DESCRIPTION  OF  WATER- VAPOR  RADIOMETERS 
Design  and  manufacturing  of  WVR's  specifically  for  use  in  geodetic  VLBI 
experiments  started  in  the  mid  seventies  at  the  Jet  Propulsion  Laboratory  (JPL)  in 
Pasadena,  CA,  where  seven  WVR’s  were  built  [Resch  el  al.  1985]  and  at  the  Onsala 
Space  Observatory  in  Sweder  where  one  WVR  was  built  [ Elgered  and  Lundh  1983]. 
Later  four  of  the  JPL  instruments  (known  as  R-series  WVR's)  were  upgraded  and 


a  new  more  compact  WVR  was  made  at  JPL  (J-series  WYR)  [Janssen  1 9 S 5] .  More 
J-series  \V\’R‘s  are  cow  being  built  and  another,  independently  designed,  WYR 
is  being  built  at  the  Geodetic  Institute  in  Bonn  for  the  YLBI  station  in  Wettzell, 

FRG. 

The  different  WVR's  that  have  been  used  in  Mark-Ill  YLBI  experiments  are 
briefly  described  in  Table  1.  System  noise-temperatures  are  all  about  600  K.  The 
WVR’s  are  all  fully  steerable  in  azimuth  and  elevation,  but  the  slew  speeds  are 
quite  different.  Even  though  all  antennas  used  for  YLBI  observations  typically 
slew  at  0.5-2° /s.  it  is  an  advantage  to  have  a  higher  slewing  speed  for  the  WYR. 
Instrumental  calibrations  are  generally  formulated  as  corrections  to  reference  loads 
or  noise-diodes  and  obtained  by  frequently  performing  elevation  scans  (also  known 
as  “tip-curves”).  This  procedure  should  be  more  successful  if  the  WYR  slews  faster 
than  the  YLBI  antenna  thus  allowing  time  between  the  YLBI  observations  for  tip- 
curves. 

The  tip-curve  method  is  sensitive  to  any  inhomogeneities  in  the  atmosphere: 
but,  provided  that  tip-curves  are  carried  out  at  diffrent  azimuth  angles,  even  sim¬ 
ple  gradients  are  easily  detected  and  modeled.  If  the  atmosphere  is  very  inho¬ 
mogeneous,  which  is  often  the  case  when  significant  amounts  of  liquid  water  are 
present,  the  noise  in  the  tip-curve  data  becomes  very  apparent  and  the  data  can 
be  downweighted  before  using  them  in  the  calibration  procedure. 
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WET  DELAYS  INTERRED  FROM  WVR  DATA 
AND  FROM  THE  KALMAN  FILTER 
When  the  Kalman  filter  is  used  to  estimate  the  atmospheric  delay  [ Herring 
et  al,  1988],  the  normal  procedure  is  to  use  the  measured  ground  pressure  to 
calculate  the  hydrostatic  delay  using  (1)  and  (2),  and  to  estimate  an  additional 
delay  which  is  then  assumed  to  be  equal  to  the  wet  delay.  The  estimated  wet  delay- 
using  the  Kalman-filter  will,  therefore,  have  an  additional  uncertainty  arising  from 
errors  in  the  a  priori  estimates  of  the  hydrostatic  delay.  Since  the  water  vapor  has 
a  different  distribution  with  height  than  has  the  dry  air,  a  special  “wet  mapping 
function”  has  to  be  used  to  calculate  the  partial  derivatives  in  the  estimation 
process.  We  have  used  the  mapping  function  presented  by  Chao  [1972]  for  the 
elevation  dependence  of  the  wet  delay. 

When  comparing  the  two  methods,  one  limitation  is  the  unavailability  of 
WVR  data  at  all  sites  and  times.  Figures  2-4  show  the  equivalent  zenith  wet  delay 
inferred  from  WVR  data  and  estimated  from  the  YLBI  data  themselves,  using  the 
Kalman  filtering  technique.  The  data  are  from  different  sites  and  involve  different 
WVRs.  We  have  omitted  the  errorbars  for  the  WVR  data  which  vary  between 
5  and  8  mm  for  the  old  R-series  WVR  and  between  2  and  4  mm  for  the  other 
instruments.  The  variations  for  a  given  instrument  are  mainly  determined  by  the 
elevation  angle  of  the  observation.  In  addition  to  these  errors,  there  are  biases  in 
the  measurements  and  in  the  inversion  algorithm  which  have  to  be  included. 

In  Figure  2  the  WVR  data  are  from  the  old  R-series  used  at  the  Haystack 
observatory,  MA;  Figure  3  shows  the  same  type  of  comparison  for  the  Mojave 
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site,  located  in  the  Mojave  desert  in  California  Figure  4  shows  one  of  the  best 
agreements  obtained.  It  consists  of  three  contiguous  experiments:  a  24  hour  long 
“Cross-Atlantic”  experiment  including  Westford,  Onsala,  and  Wettzell;  a  24  hour 
long  “IRIS”one,  including  Westford,  Fort  Davis,  Richmond,  Onsala,  and  Wettzell; 
and  a  30  hour  long  “Polar”  experiment  including  Kashima,  Fairbanks,  Mojave, 
Westford,  Onsala,  and  Wettzell.  The  Cross-Atlantic  and  Polar  experiments  are 
planned  and  scheduled  within  NASA’s  Crustal  Dynamics  Project  and  the  IRIS 
experiment  by  the  National  Geodetic  Survey. 

It  is  of  course  possible  to  modify  the  parameters  defining  the  Markov  pro¬ 
cess  in  order  to  reflect  that  the  expected  wet  delay  variations  vary  with  season 
and,  especially,  that  there  are  large  climatological  differences  between  YLBI  sites 
[ Herring  et  al.  198S]. 

ON  THE  ABSOLUTE  ACCURACY  OF  THE  TOTAL  DELAY 
The  comparisons  shown  in  Figure;  2-4  often  exhibit  similar  short-term  vari¬ 
ations  for  the  wet  delay,  but  with  a  long-term  bias  evident.  We  have  studied 
these  biases  for  the  Onsala  site  using  a  set  of  101  Mark-Ill  \  LBI  experiments  in 
which  WVR  and  VLBI  data  are  both  available  at  the  Onsala  site  for  more  than 
approximately  half  of  each  experiment.  These  experiments  were  carried  out  from 
July  1980  to  June  1987.  The  WVR  data  were  used  with  the  algorithm  presented 
by  Johansson  ei  al  [1987]  to  estimate  the  wet  delay.  We  used  ground  pressure 
measurements  together  with  the  WVR  data  to  determine  the  a  priori  atmospheric 
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delay  at  ODsala,  but  allowed  for  an  estimate  of  a  constant  correction  to  the  equiv¬ 
alent  delay  in  the  zenith  direction  for  each  experiment.  At  all  other  sites  the 
atmospheric  delays  were  estimated  by  assuming  the  random  behaviour  could  be 
represented  by  an  integrated  random  walk  process.  The  docks  at  all  sites  were 
estimated  assuming  a  random  walk.  The  statistical  parameters  of  the  Markov  pro¬ 
cess  for  the  atmospheric  delay  estimation  were  obtained  from  previous  analyses  of 
the  delay-rate  residuals  for  each  experiment  [Herrtng  ei  al ,  1988].  Ideally  each 
estimate  of  this  additional  delay  in  the  zenith  direction  would  be  zero.  There  are 
several  sources  of  error  which  will  influence  the  result: 

1.  Error  in  the  inversion  algorithm  used  with  the  WYR  data,  due  to  approxima¬ 
tions  of  the  atmospheric  profiles  of  pressure,  temperature,  and  humidity  (2  mm 
RMS  in  the  zenith  direction  for  the  algorithm  used  in  this  data  set). 

2.  Error  in  the  inversion  algorithm  used  with  the  WYR  data  due  to  uncertainties 
in  the  attenuation  coefficients  of  water  vapor  (approximately  4— 89c  of  the  wet 
delay). 

3.  WYR  instrumental  error  (RMS  scatter  3  mm.  bias  perhaps  up  to  1  cm). 

4.  Uncertainty  of  the  wet  refractivity  (the  constants  in  eq.  (3),  RMS  17c  of  the 
wet  delay). 

5.  Error  in  the  total  pressure  measurement  (estimated  peak-peak  error  at  Onsala: 
±  1  mbar,  corresponding  to  ±  2.3  mm  in  equivalent  zenith  delay). 

6.  Uncertainty  in  equation  (i)  for  the  hydrostatic  delay  (RMS  0.17c,  corresponding 
to  2.3  mm  in  equivalent  zenith  delay). 
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7.  Violation  of  hydrostatic  equilbrium  in  the  atmosphere.  Such  errors  should, 
however,  only  be  important  (larger  than  1  mm  in  equivalent  zenith  delay)  when 
strong  winds  exist.  In  this  data  set,  there  are  no  data  taken  at  Onsala  before 
March  1987  when  the  wind  has  exceeded  13  m/s  due  to  a  wind-speed  limit  for 
antenna  operation. 

8.  Any  unmodeled  effect  in  the  VLBI  data  which  is  absorbed  in  the  atmospheric 
delay  estimation  such  as  errors  in  the  mapping  function  used  to  estimate  equivalent 
zenith  delays;  these  are  mainly  sensitive  to  observations  made  at  low  elevation 
angles. 

Three  sets  of  solutions  are  presented  in  Figure  5.  Each  solution  uses  a  dif¬ 
ferent  elevation  cut-off  angle  for  the  observations  made  at  Onsala.  The  effect  of 
an  error  in  the  mapping  function  is  expected  to  be  larger  at  low  elevations,  which 
should  be  reflected  in  the  estimated  mean  bias.  Of  course  the  uncertainty  of  es¬ 
timated  zenith  bias  increases  rapidly  with  increasing  cut-off  angle  but  it  is  not 
possible  to  explain  the  bias  mainly  with  errors  in  the  mapping-function. 

Since  most  of  the  errors  associated  with  the  wet  delay  are  relative  errors 
they  can  not  alone  explain  the  results  in  Figure  5.  Such  errors  would  increase  in 
the  summer  when  the  wet  delay  is  large  and  decrease  in  the  winter.  However,  an 
instrumental  bias  of  the  WVR  could,  at  least  partly,  be  responsible  for  the  mean 
error  of  about  1  cm.  The  existence  of  an  instrumental  bias  can  be  investigated  by 
a  side-by-side  comparison  using  another  WVR  and  such  an  experiment  is  planned 
at  Onsala  in  the  near  future. 
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Another  source  of  error  that  deserves  comment  is  the  uncertainty  in  the 
observed  total  ground  pressure.  The  observed  pressures  at  the  Onsala  site  used  in 
this  data  set  have  been  compared  to  the  results  from  other  pressure  sensors  in  the 
area  and  corrected,  if  necessary.  Thereafter,  the  pressure  was  referred  to  the  height 
above  sea-level  where  the  signals  are  received  by  the  radio  telescope.  We  believe 
this  procedure  has  resulted  in  an  accuracy  of  the  ground-pressure  measurements 
of  about  1  mbar  corresponding  to  an  equivalent  zenith  delay  uncertainty  of  about 
2  mm.  Therefore,  an  effect  due  to  measurement  errors  of  the  ground  pressure  can 
not  alone  be  responsible  for  the  biases  shown  in  Figure  5. 


USING  WVR  DATA  AND  LOW  ELEVATION  ANGLE  OBSERVATION 
Low  elevation  angle  observations  are  not  important  when  delay  corrections 
are  not  estimated  ( Herring  1986).  Such  observations  will  actually  degrade  the 
accuracy  of  the  estimated  parameters  if  a  bias  error  is  present  in  the  a  priori 
delay. 

Independent  of  the  source  of  the  bias  discussed  in  the  previous  section,  it 
is  important  to  study  the  repeatability  of  estimated  baseline  lengths  for  different 
elevation  cut-off  angles  when  no  atmospheric  delays  are  estimated  at  the  sites  with 
WVRs.  Some  results  are  shown  in  Figure  6  where  again  the  set  of  101  experiments 
involving  WVR  data  at  Onsala  are  used.  For  each  set  of  solutions  made,  for  a 
given  elevation  cut-off  angle,  the  weighted  root-mean-square  (WRMS)  residuals  of 
baseline  lengths  about  the  estimated  slopes  are  presented.  The  WRMS  value  is 
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used  as  a  measure  of  repeatability.  There  is  a  small  improvement  in  repeat abiliu 


■ 
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for  all  baselines  (excluding  Onsala-Haystack)  when  we  discard  low  elevation  data. 
The  effect  is  larger  for  the  long  baselines  since  the  error  made  in  the  local  vertical 
coordinate  affects  the  baseline  length  more  in  these  cases.  The  optimum  cut-off 
angle  is  not  identical  for  the  different  baselines.  Since  it  is  an  advantage  not  to 
delete  more  data  than  necessary,  a  cut-off  angle  of  20s  seems  reasonable  in  all 
cases. 

In  order  to  show  the  effect  more  clearly,  we  have  used  the  ground  based 
model  instead  of  the  WVR  data;  the  results  are  presented  in  Figure  7.  The 
expected  accuracy  of  this  model,  averaged  over  a  year,  is  about  2  cm  RMS  in 
equivalent  zenith  delay  for  the  Swedish  west  coast  climate.  In  this  case — since 
the  expected  bias  during  each  experiment  is  larger — an  elevation  cut-off  angle  of 
about  30°  is  recommended.  Note  that  the  data  deleted  in  this  study  are  only  those 
for  observations  made  below  the  elevation  angle  cut-off  at  Onsala.  If  the  cut-off 
criteria  had  been  applied  at  other  sites  as  well  we  would  have  had  a  much  larger 
data  loss  (especially  for  the  long  baselines),  implying  a  more  rapid  increase  of  the 
WRMS  for  higher  cut-off  angles. 


RESULTS  OF  ESTIMATED  BASELINE  LENGTH  REPEATABILITIES 
The  set  of  101  experiments  was  analyzed  several  more  times,  each  time  with 
a  different  method  used  to  correct  for  the  wet  delay  at  Onsala,  but  with  the 
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atmospheric  delays  for  all  other  sites  (as  well  as  the  clocks  for  all  sites)  modeled 
as  Markov  processes. 

The  estimated  slopes  of  the  baselines  from  these  analyzes  are  presented  in 
Table  2  as  all  the  WRMS  scatters  about  these  slopes.  Based  on  the  results  pre¬ 
sented  in  Figures  6  and  7  an  elevation  cut-off  angle  was  applied  when  either  WVR 
data  or  the  ground-based  model  was  used,  and  no  estimates  were  made  for  the 
atmospheric  delay  at  Onsala.  The  obtained  repeatability  is  about  the  same  in  all 
cases  which  involve  either  the  use  of  WVR  data  as  a  priori  information  or  the 
estimation  of  the  atmospheric  delays  based  on  their  being  Markov  processes.  The 
estimated  rates  of  change  of  the  baseline  lengths  are  internally  consistent,  given 
their  uncertainties,  for  all  methods  excluding  that  using  the  ground  based  model 
with  no  further  estimation  being  made  of  the  delay  at  Onsala.  For  two  base¬ 
lines  (Onsala-Haystack  and  Onsala-Owens  Valley)  it  appears  better  to  estimate 
a  correction  to  the  zenith  delay  inferred  from  the  WVR  data  than  to  discard  low 
elevaticn  observations.  In  the  case  of  Haystack  we  know  from  Figure  6  that  the 
low  elevation  data  actually  yielded  a  lower  WRMS,  even  with  no  atmospheric  bias 
being  estimated.  Note  that  the  90%  confidence  intervals  are  relatively  large  for 
these  two  baselines  and  that  they  have  relatively  more  old  da’a  compared  to  the 
other  baselines  (see  the  mean  time  epochs  given  in  Table  2).  The  WVR  data,  as 
well  as  the  VLB1  data,  are  believed  to  be  least  accurate  in  1980. 

The  change  of  accuracy  with  time  is  further  illustrated  in  Figure  8  where 
we  have  combined  all  the  results  from  Onsala  to  Haystack  and  to  Westford.  The 


WRMS  for  the  last  18  experiments  is  8  mm,  when  the  WYR  data  are  used,  and 
11  mm  when  instead,  the  delays  are  estimated. 


CONCLUSIONS 

The  use  of  WVR  data  and  the  Kalman-filter  for  the  “wet”  delay  calibration 
give  comparable  repeatability  of  estimated  baseline  lengths  when  used  for  the 
Onsala  VLBI  site.  We  find  from  our  data  set  that  the  “best”  elevation  cut-off 
angle  when  WYR  data  are  used  is  about  20".  Low  elevation  angles  are,  however, 
important  since  they  imply  a  better  geometry  and  are  necessary  to  obtain  accurate 
results  if  atmospheric  delays  are  to  be  estimated  [Dotir  et  al  1988],  Such  estimates 
may  be  necessary  if.  for  example,  a  WYR  is  not  available  at  all  sites  or  the  WYRs 
may  not  be  working  properly  at  all  sites.  E'en  if  a  WYR  is  working  properly,  it 
may  rain  during  a  partior.  of  an  experiment,  thus  making  the  WYR  data  useless  for 
our  application.  It  is.  therefore,  of  the  greatest  importance  to  minimize  possible 
biases  in  the  total  atmospheric  delay  inferred  from  ground  pressure  and  WYR  data 
so  that  lower  elevation  angles  can  be  used  without  degradation  of  the  accuracy. 

Most  of  the  results  in  this  paper  involve  WYR  data  taken  with  one  specific 
instrument  operating  in  the  specific  climate  of  the  Swedish  west  coast.  This  WVR 
is  not  a  state  of  the  art  instrument  and  the  variations  in  the  wet  delay  at  Onsala 
are  expected  to  be  smaller  than  those  present  at  many  other  geodetic  YLBI  sites. 
Therefore,  a  more  accurate  instrument  at  a  more  humid  site  should  imply  larger 
improvements  when  a  WYR  is  used  to  correct  for  the  wet  delay. 
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Given  the  repeatability  obtained  using  the  Kalman-filter  technique  and  the 
cost  of  a  WVR,  it  may  prove  useful  to  have  a  WVR  only  at  sites  frequently  used 
and/or  where  the  expected  wet  delay  variations  are  large  The  WVR  data  can 
then  also  be  used  to  check  simultaneous  Kalman  filter  estimates  of  the  wet  delav. 
to  guard  against  other  unmodeled  errors  in  the  YLBI  data  being  absorbed  into 
atmospheric  delay  estimates:  The  WVR  is,  over  all,  an  excellent  instrument  for 
supplying  continuous,  accurate,  and  independent  information  on  the  wet  delay  in 
the  atmosphere. 
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TABLE  1.  WATER-  WAP  OR  RADIOMETERS 
USED  IN  THE  MARK-III  VLBl  EXPERIMENTS 


R-series1 

New  R-series1 

J-series1 

ASTRID2 

Frequencies 

(GHz] 

20.7  31.4 

20.7  31.4 

20.7  22.2  31.4 

21.0  31.4 

Antenna 
beam-width  [°] 

t  l 

1  l 

9  9  7 

6  6 

Reference  load 
temperatures  [K] 

313/413 

313/413 

313  + 
noise-diode 

77/313  or 
313/3603 

IF  Bandwidth 
(DSB)  [MHz] 

5-100 

10-110 

40-200 

5-500 

Slewing  speed 
AZEL  [7s] 

1.5.  1.5 

7.5.  6.6 

12.  60 

1.7.  1.7 

1  See  text 

2  The  WVR  at  the  Onsala  site  (Atmospheric  Sky  Temperature  Radiometer  for 


Interferometric  Delay  corrections) 

3  Referred  to  as  “cold"  or  “hot’'  mode.  The  WYR  was  running  in  hot  mode  during 


TABLE  2.  BASELINE  LENGTH  REP  EAT  ABILITY 


Baseline 

Onsala 

to 

Method  used  to  correct  for 
the  wet  delay  at  Onsala1 

A  priori  Adjustment 

Slope2 

(mm/year) 

WRMS3 

about  slope 
(mm) 

Wettzell 

None 

Bias 

-2.9m 

0.9 

5.7 

(55  experiments, 

None 

Markov 

-2.7± 

0.6 

4.2 

919  km,  mean 

Model(30°) 

None 

-2.2± 

0.9 

5.8 

epoch  1986.3) 

Model 

Bias 

-2.6± 

OS 

5.6 

Model 

Markov 

-2.6± 

0.7 

4.5 

WVR  (20°) 

None 

-1.9± 

0.6 

4.2 

WVR 

Bias 

-2.3± 

0.7 

4  4 

WVR 

Markov 

-2.1± 

0.7 

4.4 

Haystack 

None 

Bias 

14. 5± 

1.3 

14.1 

(30  experiments, 

None 

Markov 

15. 2± 

1.2 

12.1 

5600  km,  mean 

Model(30°) 

None 

14.3m 

2.7 

26.4 

epoch  1984.0) 

Model 

Bias 

12  0m 

1.8 

18.5 

Model 

Markov 

15. 4± 

12 

11.3 

WVR  (20°) 

None 

16. 3± 

1.3 

13.3 

WVR 

Bias 

14. 6m 

1.2 

11.5 

WVR 

Markov 

15.4m 

1.3 

12.0 

West  ford 

None 

Bias 

13.4m 

1.9 

22.2 

(78  experiments. 

None 

Markov 

13  4m 

13 

14.2 

5601  km,  mean 

Model(30c) 

None 

15. Cm 

2  1 

24.0 

epoch  1956.2) 

Mode] 

Bias 

12.7m 

1  8 

20.5 

Model 

Markov 

13  2m 

1  2 

14.2 

WVR  (20c) 

None 

12.7m 

1.1 

12.8 

WVR 

Bias 

13. 8± 

1.3 

14.5 

WVR 

Markov 

14.2m 

1.2 

13  5 

Richmond 

None 

Bias 

0.5m 

7  4 

32.1 

(24  experiments, 

None 

Markov 

0.2m 

5.6 

24.9 

7307  km,  mean 

Model(30° ) 

None 

14. 0i 

9.3 

39.7 

epoch  1986. 3) 

Model 

Bias 

0.2m 

6  1 

26.7 

Model 

Markov 

0.6m 

59 

26  2 

WVR  (20°) 

None 

4.3m 

5.1 

22  4 

WVR 

Bias 

6  3m 

5  7 

24.7 

WVR 

Markov 

6.4m 

5.8 

25.2 

OVRO 

None 

Bias 

16.9m 

2.9 

35.9 

(27  experiments, 

None 

Markov 

16. 1± 

2.6 

30.4 

7914  km,  mean 

Model(30c ) 

None 

10.4± 

3.4 

40.2 

epoch  1984.2) 

Model 

Bias 

15  0m 

2.4 

29  7 

Model 

Markov 

15.9m 

2.6 

30.4 

WVR  ,20°) 

None 

17. 7± 

3.0 

35.0 

WVR 

Bias 

15. 2± 

2.1 

24.1 

WVR 

Markov 

16.6m 

2.7 

29.4 
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* 

GRAS 

None 

Bias 

14.4- 

2.7 

29  6 

(64  experiments, 

None 

Markov 

12. 9± 

16 

26  3 

7911  km,  mean 

Me  iel(30°) 

None 

13.8= 

2.6 

363 

epoch  1965  4) 

Model 

Bias 

13.7= 

2.4 

394 

1 

Model 

Markov 

12.9= 

16 

26  1 

r 

WVR  (20°) 

None 

13.  S± 

1.9 

27.6 

i 

WVR 

Bias 

12.9= 

1.9 

26. 2 

WVR 

Markov 

136= 

1.3 

27.3 

1  The  o  prtori  information  is  none,  or  the  ground-based  model, or  WVR.  If  an  ele¬ 
vation  cut-off  angle  is  applied,  its  value  is  given  in  paraDtheses  If  an  atmospheric  delay  is 
estimated,  it  is  assumed  to  be  either  a  constant  in  the  zenith  direction  (bias)  or  a  Markov 
process. 

2  The  sigmas  are  scaled  so  that  reduced  x 2  Per  degree  of  freedom  =  1  ( Herring  et 
al  1966). 

3  weighted  root-mean-square 
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FIGURE  CAPTIONS 


Figure  1 

WVR  mesurements  carried  out  during  four  days  in  November  19S6.  T^e 
large  errors  in  the  inferred  wet  delays  from  \\  V'R  data  taken  during  rain  is  in  this 
case  further  increased  by  water-drops  forming  on  the  covers  of  the  horn  antennas. 
The  errorbars  for  the  WVR  data  have  been  omitted. 

Figure  2 

Simultaneous  WVR  measurements  and  Kalman  filter  estimates  of  the  equiv¬ 
alent  zenith  wet  delay  at  the  Haystack  Observatory.  The  old  R-series  WVR  is 
used. 

Figure  3 

Simultaneous  WVR  measurements  an  -1  Kalman  filter  estimates  of  the  equiv¬ 
alent  zenith  wet  delay  at  the  Mojave  VLB1  site  The  new  R-series  WVR  is  used. 

Figure  4 

Simultaneous  WVR  measurements  and  Kalman  filter  estimates  of  the  equiv¬ 
alent  zenith  wet  delay  at  the  Onsala  Space  Observatory.  The  ASTRID  WVR  is 
used. 
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Figure  5 


A  constant  atmospheric  delay  at  Onsala  estimated  in  the  addition  to  an  a 
prion  delay  consisting  of  the  hydrostatic  delay  calculated  using  the  total  ground 
pressure  and  the  wet  delay  inferred  from  WVR  data.  The  estimation  is  done 
for  101  VLB1  experiments  three  times,  each  time  with  a  different  elevation  cut-off 
angle  at  Onsala  in  order  to  check  the  sensitivity  of  the  estimated  value  for  mapping 
function  errors. 


Figure  6 

Weighted  RMS  residuals  of  estimated  baseline  lengths  about  a  best  fit  straight 
line.  The  hydrostatic  delay  and  WVR  data  are  the  a  priori  information  used  and 
no  estimate  of  the  atmospheric  delay  at  Onsala  is  done.  Eight  sets  of  solutions  are 
made,  each  one  with  a  different  elevation  cut-off  angle  at  Onsala.  The  baselines 
are  from  Onsala  to:  Wettzell  (2.  919  km,  55  experiments).  Haystack  (K,  5600  km. 
30  experiments),  Westford  (E,  5601  km,  78  experiments),  Richmond  (R.  7307  km. 
24  experiments),  OYRO  (O.  7914  km.  27  experiments),  and  GRAS  (G.  7941  km. 
64  experiments). 


Figure  7 

Weighted  RMS  residuals  of  estimated  baseline  lengths  about  a  best  fit  straight 
line.  The  hydrostatic  delay  and  the  ground  based  model  are  the  o  priori  infor¬ 
mation  used  and  no  estimate  of  the  atmospheric  delay  at  Onsala  is  done.  For  the 
explanation  of  the  baselines  see  the  caption  to  Figure  6. 
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Figure  6 


Estimated  baseline  lengths  from  Onsala  to  the  “combined’'  Haystack/West- 
ford  site.  The  circles  denote  measurements  to  Haystack  ar.d  the  triangles  denote 
measurements  to  Westford.  Either  Haystack  or  Westford  or  both  have  been  in¬ 
volved  in  all  the  101  experiments  studied.  Since  both  sites  have  been  involved  in 
several  experiments  simultaneously  there  are  108  measurements  all  together.  It 
is  seen  that  both  the  formal  error  of  the  individual  experiments  and  the  repeata¬ 
bility  have  improved  with  time.  The  following  milestones  should  be  noted:  June 
1982:  The  Mark-Ill  system  is  upgraded  from  7  to  1*3  videoconverters  implying  bet¬ 
ter  group  delay  measurements;  May  1SS5:  Westford  installs  cooled  receiver  and 
replaces  Haystack  in  almost  all  CDP-experiments.  August  1980:  Onsala  installs 
cooled  receiver,  and  March  1987:  Onsala  installs  a  dual  frequency  feed  in  the  20m 
radome  enclosed  telescope  which  previously  was  used  for  the  X-band  only. 
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VLBI  geodesy:  2  parta-per-billion  precision  in  length  determinations  for  transconti¬ 
nental  baselines 


J.L.  Davis,  T.A.  Herring,  and  I.I.  Shapiro 
Harvard-Smithsonian  Center  for  Astrophysics 
Cambridge,  Massachusetts  02138 


ABSTRACT.  We  have  used  very-long-baseline  interferometry  (VLBI)  to  make  twenty-two  indepen¬ 
dent  measurements,  between  September  1984  and  December  1986,  of  the  length  of  the  3900-km 
baseline  between  the  Mojave  site  in  California  and  the  Haystack/Westford  site  in  Massachusetts. 
These  experiments  differ  from  the  typical  geodetic  VLBI  experiments  in  that  a  large  fraction  of 
observations  are  obtained  at  elevation  angles  between  46  and  10°.  Data  from  these  low  elevation 
angles  allows  the  vertical  coordinate  of  site  position,  and  hence  the  baseline  length,  to  be  estimated 
with  greater  precision.  For  the  sixteen  experiments  processed  thus  far,  the  weighted  roct-mean- 
square  scatter  of  the  estimates  of  the  baseline  length  is  8  mm.  We  discuss  these  experiments,  the 
processing  of  the  data,  and  the  resulting  baseline  length  estimates. 


2.  INTRODUCTION 

In  recent  years,  the  precision  of  estimates  of  baseline  length  obtained  from  the  analysis  of  Mk- 
III  very-long-baseline  interferometry  (VLBI)  data  has  passed  below  the  level  of  10  parts- per-biilion. 

The  limiting  source  of  error  continues  to  be  the  atmosphere.  The  effects  of  atmospheric  errors  on 
estimates  of  baseline  length  are  discussed  in  Lanyi  [1984),  Davis  et  of.  [1985],  Davis  [1956],  Herring 
[1986],  and  Treuhaft  and  Lanyi  [1987].  The  typical  treatment  of  atmospheric  propagation  delay  in 
the  analysis  of  geodetic  VLBI  data  involves  the  use  of  surface  meteorological  measurements  (and  ( 

possibly  radiometric  measurements;  see  Elgered  tt  aJ.  [1987])  to  determine  an  a  priori  value  for 
the  zenith  propagation  delay,  and  the  use  of  a  mapping  function.  (The  mapping  function  describes 
the  elevation-angle  dependence  of  the  propagation  delay,  ».e.,  the  number  of  air  masses  traversed.) 

Corrections  to  the  a  priori  value  of  the  zenith  delay  are  then  estimated,  along  with  the  relevant 
geodetic  parameters,  from  the  VLBI  data. 

The  unique  elevation-angle  dependence  of  the  atmospheric  propagation  delay  allows  a  precise  I 

estimate  to  be  made  of  the  senith-delay  corrections.  However,  the  sensitivity  of  the  VLBI  gToup- 
delay  measurement  to  a  change  in  the  tenith  delay  is  highly  correlated  with  the  sensitivity  of  the 
group  delay  to  a  change  in  the  vertical  coordinate  of  site  position  for  small  ranges  of  air  mass. 

By  observing  sources  at  low  elevation  angles,  we  can  decrease  the  correlation  of  the  estimates  of 

the  senith-delay  corrections  with  those  of  the  vertical  coordinate  of  site  positions,  and  thereby 

obtain  a  more  precise  estimate  of  those  (and  other)  parameters.  Nevertheless,  many  schedules  for  I 

geodetic  VLBI  experiments  include  no  observations  below  about  10s  elevation  at  any  site.  This 

omission  is  primarily  intended  to  avoid  the  effects  of  errors  in  the  mapping  function  used,  since  such 

errors  are  generally  worse  at  lower  elevations.  (It  is  harder  to  predict  atmospheric  properties  at 

greater  horizontal  distances  from  the  site.)  W'e  therefore  pose  the  question:  What  is  the  “optimum' 

minimum  elevation  angle  to  minimize  errors  in  baseline-length  estimates? 
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Pigure  1.  Estimates  of  the  length  of  the  baseline  from  Haystack/VVestford  to  Mojave  (see  text)  A 
mean  value  of  3,904,144,248  mm  has  been  subtracted  from  the  values  shown. 

2.  DESCRIPTION  OF  EXPERIMENTS  and  RESULTS 

In  order  to  study  this  question,  we  designed  YLBI  experiments  known  as  “low-elevation”  ex¬ 
periments  in  reference  to  their  scheduling  strategy:  the  attempt  is  made  to  observe  low  in  the  sky  as 
frequently  as  possible,  within  the  limits  of  the  antennas  used.  For  the  standard  set  of  low. elevation 
experiments,  we  used  the  Mojave  site  (em,n  =  8°)  in  California  and  the  Haystack/ Westford  site 
(«m,n  =  4°)  in  Massachusetts. 

The  processing  of  the  data  from  these  experiments  differed  in  several  respects  from  the  pro¬ 
cessing  which  we  employed  several  years  ago  (see,  e.j.,  Clark  ef  cl.  (1985 j)  for  geodetic  VLBI  data 
A  new  mapping  function  for  the  dry  atmosphere,  developed  by  Davis  et  cl.  (1965],  and  accurate 
to  about  10  mm  at  5°  elevation  for  a  wide  range  of  atmospheric  conditions,  was  used.  Another 
improvement  was  the  use  of  a  Kalman  filter  to  estimate  clock  offsets  and  renith  delay  corrections. 
Each  of  these  corrections  is  modeled  as  a  time- varying  stochastic  process,  the  statistics  of  which  are 
deduced  from  other  sources;  the  values  of  the  corrections  for  the  epochs  of  the  VLBI  observations 
are  estimated  by  the  Kalman  filter.  Finally,  we  used  radiometric  data  to  estimate  the  “wet  delay’' 
when  these  data  are  available. 

In  Figure  1,  we  present  the  estimates  of  the  Haystack-Mojave  baseline  lengths  as  a  function  of 
experiment  date.  The  error  bars  represent  the  standard  deviation  of  the  estimate  of  the  baseline 
length,  deduced  from  the  propagation  of  the  stochastic  errors  involved.  These  errors  include  the  the 
measurement  errors  of  the  group-delay  observations  and  the  stochastic  behavior  of  the  clocks  and 
atmospheres.  The  weighted  root-mean-square  scatter  of  the  estimates  about  their  weighted  mean 
is  8  mm,  representing  a  fractional  repeatability  of  2  parts-per-billion. 

This  work  was  supported  by  Air  Force  Geophysics  Laboratory  contract  F-19628-86-K-0025: 
NASA  grant  NAG5-538,  and  NSF  grants  EAR-83-06380  and  EAR-86-189S9. 
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VLB!  studies  of  the  nutations  of  the  earth 


T.  A.  Herring.  Harvard-Smithsonian  Center  for  Astrophysics 
60  Garden  Street 
Cambridge,  MA.  02138 


ABSTRACT.  The  application  of  very-long-baseline  interferometry  (YLB1)  to  the 
study  of  the  nutations  of  the  earth  has  yielded  unprecedented  accuracy  for  the 
experimental  determination  of  the  coefficients  of  the  nutation  series.  The  anal¬ 
ysis  of  six  years  of  VLBI  data  has  yielded  corrections  to  the  coefficients  of  the 
seven  largest  terms  in  the  1AU  1980  nutation  series  with  periods  of  one  year  or 
less,  with  accuracies  approaching  the  truncation  error  of  this  nutation  series  (0.1 
mas).  The  nutation  series  coefficients  computed  from  the  VLBI  data,  and  those 
obtained  from  theoretical  considerations  (the  1AU  1980  nutation  series),  are  in  ex¬ 
cellent  agreement.  The  largest  corrections  are  to  the  coefficients  of  the  retrograde 
annual  nutation  [2.0  ±  0.1  mas),  the  prograde  semiannual  nutation  '(0.5  -  t  0.4  ; 
±  0.1  mas’,  and  the  prograde  13.7  day  nutation  [-0.4  *  0.1  mas).  (The  imaginary 
term  for  the  semiannual  nutation  represents  a  term  90c  out-of-phase  with  the  ar¬ 
guments  of  the  nutation  series.)  The  geophysical  implications  of  these  results  are 
currently  under  active  investigation.  We  discuss  the  methods  used  to  extract  the 
nutation  information  from  the  VLBI  data,  the  calculations  of  the  uncertainties  o' 
the  resultant  corrections  to  the  coefficients  of  the  nutation  series,  and  the  current 
research  into  the  nutations  of  the  earth. 


1.  INTRODUCTION 

The  determination  of  the  coefficients  of  the  nutation  series  using  ver\ -long- 
baseline  interferometry  data  has  been  discussed  recently  in  a  number  of  papers 
(Herring  it  al.,  1983;  Guinn  et  al .,  1984;  Herring  ei  a/.,  1985;  Eubanks  et  a !.,  1955. 
Guinn  et  at.,  1986;  Herring  et  al.,  1986a;  Herring  el  al.,  1986b;  and  Htmwtch  et  a!.. 
1986).  A  detailed  description  of  the  analysis  techniques  used  for  determining  the 
corrections  to  the  coefficients  of  the  nutation  series  is  given  in  Herring  et  al.,  1956a 
The  interpretation  of  the  results  given  in  Herring  et  al.,  1986a  is  discussed  in  detail 
by  Guinn  et  of.,  1986.  The  interpretation  of  these  results  has  also  been  discussed 
by  Wohr  and  Bergen  el  of.,  1986;  Yoder  and  Ivinus,  1986;  and  Yoder  and  Ivmus, 
1987.  Here  we  review  briefly  the  procedures  used  to  estimate  the  corrections  to  the 
coefficients  of  the  nutations,  and  we  present  the  results  obtained  from  the  analysis 
of  412  VLBI  experiments  obtained  during  the  6.5  year  period  between  July  1950 
and  February  1987. 
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2.  DATA  ANALYSIS  TECHNIQUE 

The  techniques  we  have  used  to  estimate  corrections  to  the  IAU  1963  nuta¬ 
tion  series  have  been  discussed  in  Htrring  et  al  (1985)  and  Herring  et  al.( 1966a, 
For  each  VLBI  observing  session  (typically  of  24  hours  duration),  we  estimate 
corrections,  6 and  iAt,  to  the  nutation  angles  computed  from  the  IAU  I960 
nutation  series.  In  a  post-processing  operation,  we  use  these  nutation-angle  cor¬ 
rections  to  estimate  corrections  to  the  coefficients  of  certain  terms  in  the  nutation 
series.  An  alternative  technique  is  to  directly  estimate  the  corrections  to  coeffi¬ 
cients  of  the  nutation  series  from  the  VLB!  data  themselves.  Results  obtained 
using  this  technique  have  been  presented  by  Himvich  et  a!.,  1986. 

The  nutation-angle  corrections  are  not  only  used  to  estimate  the  corrections 
to  the  coefficients  of  the  nutation  series,  but  they  are  also  be  used  to  study  the 
noise  characteristics  of  the  nutation  estimates  (see  Herring  et  al. ,1986a:.  and  to 
study  phenomena  not  explicitly  included  in  the  standard  nutation  series  such  as 
the  “free  core-nutation"  (see  Herring  et  al .,  1986b).  In  particular,  we  will  discuss 
the  power-spectra!  density  (PSD)  function  of  the  nutation-angle  corrections.  The 
calculation  of  the  PSD  function  from  the  nutation-angle  corrections  have  beer, 
discussed  in  Htrring  et  al.,  1965.  The  technique  we  adopt  is  to  “lightly"  smooth 
the  nutation-angle  corrections  using  a  Gaussian  filter  with  a  full-width-at-half- 
maximum  (FWHM)  of  5  days,  and  to  then  compute  the  PSD  function  from,  these 
uniformly  spaced  smoothed  values  (5  day  spacing).  The  (assumed  uncorrelatec; 
noise  in  the  nutation-angle  corrections  is  propagated  through  the  smoothing  and 
the  fast  fourier  transform  (FFT)  algorithms  to  compute  the  99.591  confidence  limit 
for  the  PSD  function. 

3.  RESULTS 

The  VLBI  data  we  have  analyzed  are  from  412  observing  sessions  carried 
out  by  the  NASA  Crustal  Dynamics  Project  and  the  NGS  IRIS  program,  between 
July  1980  and  February  196T.  This  data  set  is  an  extension  of  those  discussed 
in  Herring  et  a/. (1965)  and  Herring  et  a!. (1966a  and  b).  In  Table  1,  v.e  give  the 
estimates  of  the  amplitude^  of  the  circular  nutations  that  can  be  resolved  w  ith  the 
limited  temporal  range  of  the  VLBI  data.  The  values  for  and  sir.  c.  we 
used  to  obtain  these  amplitude  estimates  are  shown  in  Figure  1.  The  PSD  function 
computed  from  the  nutation-angle  corrections  shown  in  Figure  1  is  shown  in  Figure 
2.  The  three  largest  peaks  in  the  PSD  function  correspond  to  the  three  largest 
corrections  to  the  coefficients  of  the  nutation  series  given  in  Table  1 .  There  is  also  a 
large  peak  corresponding  to  the  retrograde  frequency  closest  to  1  cycle  per  sidereal 
day.  Any  corrections  to  the  long  period  nutations  (18.6  and  9.3  years)  would  ail 
appear  in  this  term.  With  only  6  5  years  of  data  we  can  not  yet  reliably  estimate 
corrections  to  the  nutation  series  terms  with  these  long  periods.  The  remainder 
of  the  PSD  function  appears  to  be  consistent  'vjth  the  while-noise  statistics  used 
to  compute  the  99.5%  confidence  interval. 

The  results  given  in  Table  1  are  very  consistent  with  previously  published 
results  within  the  uncertainties  of  these  earlier  analyses.  The  correction  to  the 
retrograde  ennua!  nutation  has  been  interpreted  as  being  due  to  the  effects  of  the 
earth’s  fluic.  core  being  more  strongly  coupled  to  the  mantle  than  was  inferred 
for  the  calculation  of  the  IAU  1980  nutation  series  (Guinn  et  al.,  1966).  The 
corrections  to  the  prograde  semiannual  term  are  partly  due  the  coupling  of  the 
fluid  core  to  the  nmntle  and  partly  due  to  the  effects  of  ocean  tides  on  the  nutations 
( Satao  and  H’ohr,  1961).  The  correction  to  the  13.7  day  prograde  term  has  not 
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TABLE  1.  Estimates  of  Corrections  to  Amplitudes  of  Circular  Nutations  of  lAV 
1980  Nutation  Series  (  Wohr,  1961). 


Corrections*’ 


Period  Observed  In-Phase  Out-of-Pnase 

in  inertial  Amplitude,  6ar,  6a, , 

•pace,  days  mas  mas  mas 


162  6  prograde 

-548  62 

0  45 

-0  41 

162  6  retrograde 

-24  63 

-0.11 

-0 .07 

13  7  prograde 

-94  45 

-0  37 

-0  00 

13  7  retrograde 

-3  63 

-0  03 

0  05 

365  3  prograde 

25  71 

0  05 

0  16 

365  3  retrograde 

-33  13 

-2  07 

0  36 

2*  6  prograde 

14  50 

-0  02 

0  Cf 

27.6  retrograde 

-13  79 

0  03 

-0  Cl 

121.7  prograde 

-21  56 

-0  05 

-0  03 

121  7  retrograde 

-0  95 

-0  04 

-0  03 

9  1  prograde 

-12  53 

-0  10 

0  00 

9  1  retrograde 

-0  43 

-0  01 

-o  o: 

31  8  prograde 

3  19 

-0  01 

-0  03 

31  8  retrograde 

-3.02 

o  o; 

0  00 

433  2  prograde 

0.06* 

-0  06 

-0  00 

433.2  retrograde  FCN 

0  33* 

0.32 

-0  04 

‘The  standard  deviation  of  each  correction  is  estimated  to  be  0.10  mas  The 
standard  deviations  are  obtained  using  the  techniques  described  in  Htrr mg  et 
ai  (1986),  with  the  root-sum-squares  addition  of  a  contribution  of  0.0?  mas  to 
account  for  (i)  the  effects  of  the  truncation  of  the  IAU  1960  nutation  series  to  the 
nearest  0.1  mas,  (ii)  the  possible  effects  of  atmospheric  excitation  of  the  nutation, 
and  (iii)  the  effects  of  ocean  and  earth  tide  modeling  errors  in  the  analysis  of  the 
VLBI  data.  See  also  Herring  et  at.  (1986)  for  definitions  of  the  quantities  listed 
tThe  amplitudes  of  these  terms  are  the  root-sum-squares  of  the  amplitudes 
of  the  real  and  imaginary  components. 


been  satisfactorily  explained.  Since  this  nutation  term  is  the  one  most  affected  by 
the  elasticity  of  the  earth,  the  correction  may  arise  from  the  treatment  of  elasticity 
or  the  neglect  of  the  effects  of  anelasticity  in  the  IAU  nutation  series.  The  root- 
mean-square  scatter  of  the  remaining  25  corrections  given  in  Table  1  is  0.06  mas. 
which  is  consistent  with  our  estimated  uncertainties  of  0.10  mas. 

4.  CONCLUSION 

The  quality  of  the  VLBI  estimates  of  the  corrections  to  the  nutation  series 
will  improve  as  further  VLBI  experiments  are  carried  out.  The  VLBI  data  from 
the  international  earth  rotation  service  will  continue  to  provide  high  quality  mea¬ 
surements  of  the  nutation  angles  on  a  regular  basis  (once  every  5  days).  These 
data  and  those  obtained  by  the  NASA  Crustal  Dynamics  Project  will  ensure  that 
the  coefficients  of  selected  terms  in  the  nutation  series  will  be  determined  with 
precisions  of  less  than  0.1  mas.  However,  the  current  nutation  series  is  truncated 
at  0.1  mas,  and  thus,  will  soon  be  the  limiting  error  source  in  the  interpretation 
of  the  corrections  to  the  nutation  series  Since  all  modern  nutation  series  are  ob¬ 
tained  by  convolving  the  response  function  of  a  ‘‘realistic"  model  earth  with  the 
nutations  of  a  rigid  earth,  we  will  soon  start  developing  a  rigid-earth  nutation 
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Figure  1.  Estimates  of  £A<  and  /Avsintc  from  the  VLBI  data  analyzed  :r.  this 
paper.  The  solid  line  in  each  figure  is  computed  from  the  corrections  giver,  ir. 
Table  1  and  from  long-period  terms  (periods  >  IS. 6  years)  which  are  not  giver, 
in  the  table  because  their  values  are  highly  correlated  and  the  estimates  a^e  not 
considered  reliable. 
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series  with  an  accuracy  of  0  01  mas.  This  ne«  nutation  series  will  provide  ar. 
accurate  framework  in  which  the  observed  nutations  of  the  earth  can  be  compared 
with  that  computed  from  geophysical  theory 
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day  (cpsd)  is  due  to  use  of  a  (5  day  FWHM)  Gaussian  filter.  The  frequency  spa. 
ing  in  the  PSD  function  is  1/1830  cpsu.  At  the  top  of  the  figure  we  have  shov 
the  frequency  positions  of  some  of  the  major  terms  in  the  nutation  series  V 
have  also  marked  the  frequency  components  in  the  PSD  function  closest  to  the 
frequencies  with  triangles. 
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ABSTRACT.  We  analyzed  six  years  of  very-long-baseline  interferometry  (VLB!) 
data  and  determined  corrections  to  the  coefficients  of  the  seven  terms  with  the 
largest  amplitudes  in  the  IAU  1980  nutation  series.  Our  analysis  yields  results 
consistent  with  earlier  analyses  of  smaller  sets  of  VLBI  data,  within  tne  uncertain¬ 
ties  of  the  latter.  Here,  we  restrict  discussion  to  the  freely  excited  core-nutation 
or  “free  core-nutation"  (FCN).  Our  analysis  yields  an  estimate  of  0.33  ±  0.12  mas 
for  an  assumed  constant  amplitude  of  the  FCN.  which  allows  us  to  place  an  upper 
bound  on  it  of  0.6  mas  (99.5%  confidence  limit).  We  also  studied  possible  tempo¬ 
ral  variations  of  the  complex  amplitude  of  the  FCN  by  modeling  it  as  a  stochastic 
process  with  a  white  noise  excitation.  We  detected  no  statistically  significant 
variations  of  this  amplitude  for  the  six-year  interval  spanned  by  the  VLB]  data. 
However,  in  the  neighborhood  of  one  cycle  per  day,  the  power  spectral  density 
of  the  atmospheric  surface  loading  is  estimated  from  global  weather  data  to  be 
0.24  (g  cm~2)2  day,  about  five  times  larger  than  the  largest  such  power  spectra! 
density  that  would  be  consistent  with  the  upper  bound  on  the  amplitude  of  the 
FCN  placed  by  the  VLBI  data.  Thus,  we  conclude  that  this  estimate  is  too  high 
and  that,  if  the  FCN  were  excited  by  surface  loads  with  frequencies  near  one  cycle 
per  day,  then  the  power  spectral  density  of  these  loads  must  be  <0.06  (g  cm~2)! 
day  (99.9%  confidence  limit). 


1.  INTRODUCTION 

The  recent  investigations  of  the  nutations  of  the  earth  using  VLBI  data  ( Her¬ 
ring  et  d.,  1983;  Gtn'nrc  et  d.,  1984;  Herring  et  d.,  1985;  Eubanks  el  d.,  1985;  and 
Herring  et  d.,  1986)  have  all  disclosed  corrections  to  the  coefficients  of  the  terms  in 
the  IAU  1980  nutation  series  corresponding  to  the  retrograde  annual  nutation  and 
the  prograde  semiannual  nutation.  These  corrections  can  be  partially  explained  by 
changing  the  resonance  frequency  (in  a  frame  rotating  with  the  earth)  of  the  “core 
nutation"  from  -(1+  5^5}  cycles  per  sidereal  day  (cpsd)  to  -(1+  ^ )  cpsd  {Guinn 
et  d., 1986).  Other  phenomena  contributing  significantly  to  the  corrections  might 
be  ocean  tides  (  Wahr  and  Sasoo,  1981)  and  the  anelasticity  of  the  earth’s  mantle 
(Wahr  and  Bergen,  private  communication,  1986).  These  latter  two  effects  are 
not,  however,  nearly  large  enough  by  themselves  to  explain  the  ss  2  mas  correc¬ 
tion  to  the  amplitude  of  the  retrograde  annual  nutation.  These  small  corrections 
aside,  the  ove  all  agreement  between  the  VLBI  results  and  the  IAU  1980  nutation 
series  is  evidence  for  the  existence  of  the  predicted  effect  of  the  core-nutation  res¬ 
onance  on  the  earth's  rotation:  this  resonance  contributes  as  much  as  17  mas  to 
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the  nutation  series,  in  particular  to  the  semiannual  term.  Nonetheless,  the  free 
excitation  of  this  mode  has  never  been  unambiguously  detected.  Prior  analyses 
of  VLBI  data  have  placed  only  upper  bounds  of  under  1  mas  on  the  amplitude 
of  this  free  core  nutation  (FCN)  [Herring  et  of.,  1985;  Eubanks  et  oi,,  1985;  and 
Herring  et  al.,  1986]. 

In  this  paper,  we  examine  the  FCN  in  more  detail.  In  particular,  we  inves¬ 
tigate  the  possibility  that  our  single  (“average”)  estimate  of  the  amplitude  of  the 
FCN  mav  be  significantly  less  than  the  maximum  value  reached  by  this  amplitude 
during  the  six-year  interval  spanned  by  the  VLBI  data,  due  to  large  fluctuations 
in  this  (complex]  amplitude  during  this  interval.  This  situation  could  arise  if  the 
FCN  were  heavily  damped  and  highly  excited.  An  estimated  lower  bound  on  the 
damping  time,  combined  with  an  independent  estimate  of  the  expected  amplitude 
of  the  FCN  can  be  used  to  Quantify  statistically  the  possible  fluctuations  in  this 
amplitude.  We  can  obtain  a  lower  bound  for  the  damping  time  from  the  VLBI  es¬ 
timate  of  the  amplitude  of  the  out-of-phase  component  of  the  retrograde  annua! 
nutation.  An  estimate  of  the  expected  amplitude  of  the  excitation  of  the  FCN 
can  be  obtained  from  other  geophysical  data.  In  particular,  we  take  the  expected 
amplitude  of  the  FCN  from  the  magnitude  obtained  by  Eubanks  et  al.{  1985)  of 
the  power  spectral  density,  at  nearly  diurnal  periods,  of  the  P2j  component  of  the 
atmospheric  pressure  loading  of  the  earth's  surface  (see  H'a/ir  and  Sasao  (198 lj 
for  normalization  of  spherical  harmonics). 

2.  DATA  ANALYSIS  TECHNIQUE 

The  techniques  we  have  used  to  estimate  corrections  to  the  IAU  1980  nuta¬ 
tion  series  have  been  discussed  in  Herring  et  of.(1985)  and  Herring  et  o/.(1986).  For 
each  VLBI  observing  session  (typically  of  24  hours  duration),  we  estimate  correc¬ 
tions,  6  AV'  and  6  Ac,  to  the  nutation  angles  computed  from  the  IAU  1960  nutation 
series.  In  a  post-processing  operation,  we  use  these  nutation-angle  corrections  to 
estimate  corrections  to  the  coefficients  of  certain  terms  in  the  nutation  series.  The 
analysis  used  in  this  paper  is  similar  to  our  previous  analyses  except  that  here  we 
allow  our  model  of  the  complex  amplitude  of  the  FCN  to  vary  stochastically  dur¬ 
ing  the  interval  spanned  by  the  data  set.  A  Kalman  filter  (see,  e.g.,  Licbclt,  1967) 
was  used  to  estimate  simultaneously  the  time-dependent  complex  amplitude  of 
the  FCN,  and  the  corrections  to  the  coefficients  of  selected  terms  in  the  IAU  I960 
nutation  aeries. 

The  stochastic  model  used  to  represent  the  complex  amplitude  of  the  FCN 

was 


f(f  +  At)  =  f(f)e"oA*  +  $f(0--  (1) 

where,  respectively,  ff t)  and  f(t  +  At)  are  the  complex  values  of  the  FCN  at  time; 
t  and  t  +  At;  a  (>0  )  is  the  inverse  damping  time  for  the  FCN  and  stems  from 
the  imaginary,  or  dissipative,  part  of  the  FCN  resonance  frequency;  and  6g(t)  is 
the  complex  excitation,  integrated  over  the  interval  t  to  t  +  At  (At  C  a-1).  The 
corresponding  contributions  of  the  FCN  to  the  nutation  angles  are 

$A<rCtv(t)  ae  -fr(t)cOs((n  +  Wfc*)*)  “  f»(t)  *ir  ((n  +  UpcN)t)  (2)  1 

and 

$AV>rc;v(f)*vnfo  =  -f,(t)sin((n  +  w/-CN)t)  +  f.(<)  cos((D  +  wrcN)t),  (3) 

where,  respectively,  frft)  and  f,(t)  are  the  real  and  imaginary'  parts  of  J(t  1  = 

(r (t)  -«f,(t);  t  is  sidereal  time  (measured  for  convenience  from  the  epoch  J2000J;  H 
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is  the  earth’s  rotation  rate;  and  u ’fcn  i*  the  real  component  of  the  FCN  frequency. 
To  apply  a  Kalman  filter  to  the  estimates  of  6Ac(t)  and  £A v'(l)  obtained  from 
the  VLBI  data,  we  need  to  establish  the  statistical  properties  of  the  excitation  of 
the  FCN.  However,  we  do  not  know  the  excitation  mechanism  for  the  FCN.  The 
most  efficient  mechanism  seems  to  be  through  surface  loading  (Saaao  and  W  ahr, 
1981).  But,  very  little  is  known  about  the  properties  of  non-tidal  surface  loading 
with  periods  near  one  day.  Likely  candidates  for  the  excitation  are  loads  from 
atmospheric  pressure  variations,  and,  possibly,  from  sea  level  changes  (Sasoo  and 
IFohr,  1981).  Since  we  know  so  little  about  the  spectra  of  these  excitations,  we  will 
make  the  simple  assumption  that  all  of  them  have  white-noise  spectra  for  periods 
near  one  day.  We  will  characterize  the  amplitude  of  the  white-noise  excitation 
by  its  variance,  o*  =  (5f(l)£f*(l)),  where  *  denotes  complex  conjugation,  and  {/ 
denotes  expectation. 

The  magnitude  of  ff*  can  be  established  by  balancing  the  power  input  into 
the  FCN  by  the  excitation  with  that  dissipated.  From  equation  (1),  we  have  for 
steady-state  conditions  (assuming  aAf  Cl): 

=  2(fDoAf.  (4) 

Therefore,  if  we  know  a  and  the  expectation  of  the  square  of  the  magnitude 
of  the  (integrated)  complex  excitation  of  the  FCN,  we  can  compute  the  variance 
of  the  white-noise  excitation  to  be  used  in  the  Kalman  filter. 

3.  RESULTS 

The  VLBI  data  we  have  analyzed  are  from  370  observing  sessions  carried  out 
by  the  NASA  Crustal  Dynamics  Project  and  the  NGS  IRIS  program  between  July 
1980  and  August  1986.  This  data  set  is  an  extension  of  those  discussed  in  Herring 
tt  of. (1985)  and  Herring  et  al. (1986).  In  Table  1,  we  give  the  estimates  of  the 
amplitudes  of  the  circular  nutations  that  can  be  resolved  with  the  limited  temporal 
range  of  the  VLBI  data.  The  values  for  6  Ac  and  6  At£  sin  1 e  we  used  to  obtain  these 
amplitude  estimates  are  shown  in  Figure  1.  These  amplitude  results  are  based 
on  the  assumption  that  the  amplitude  of  the  FCN  was  constant  for  the  six  years 
spanned  by  the  data  set,  i.e.  that  o*~ 0,  whence  the  Kalman  filter  solution  reduces 
to  conventional  weighted-least  squares.  In  this  solution,  we  not  only  estimated 
the  amplitude  of  the  FCN,  but  also  the  amplitude  of  a  prograde  (“control')  term 
with  the  frequency  -(1  -  jyj)  cpsd.  Since  no  resonance  is  thought  to  exist  for 
any  such  prograde  term,  any  signal  at  this  frequency  is  probably  due  to  noise  or 
to  model  deficiencies,  and  could  be  used  as  a  measure  of  the  uncertainty  affecting 
the  estimate  of  the  FCN  amplitude.  The  resultant  estimates  of  the  amplitudes  of 
the  FCN  and  the  prograde  terms  are,  respectively,  0.33  ±0.12  and  0.04  i  0.1 2  mas. 

The  amplitude  of  the  FCN  can  also  be  estimated  from  other  data.  Eubanks 
ef  al.  (1985)  have  estimated  from  global  weather  data  that  the  power  spectral- 
density  (PSD)  of  the  Pji  component  of  the  atmospheric  pressure  field  is  0.24  (gm 
cm-J)J  day.  (No  uncertainty  for  this  value  is  given  by  the  authors.  However,  the 
estimate  is  stated  to  be  “conservative.")  This  PSD, combined  with  the  lower  bound 
for  damping  time  deduced  from  the  V'LBI  data  (see  below),  suggests  that  the 
FCN  amplitude  should  be  >0.8  mas  (see  Sasao  and  W'ahr,  1981  for  computation 
methods).  From  this  lower  bound  for  the  FCN  amplitude,  and  the  upper  bound  for 
a  obtained  from  the  out-of-phase  component  of  the  retrograde  annual  nutation, 
we  estimate  -  lower  bound  for  o* .  Taking  a  to  be  <  2*  x  3.3x10“*  d-1  (99.5CT 
confidence  limit)  [gee  Gwnn  el  al. (1986)  for  computation  methodj,  we  find  c *  > 
5.8  x  10-4At  mas',  where  Al  is  in  days.  This  lower  bound  for  o *  indicates 
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TABLE  1.  Estimates  of  Corrections  to  Amplitudes  of  Circular  Nutations  of  1AU 
1980  Nutation  Series  (WoAr,  1981). 

Corrections* 


Period 
in  inertial 
space,  days 

Observed 

Amplitude, 

KQtl 

In-Phase 

oaT, 

PM 

Out-of-Phase 

«a„ 

m&B 

182  6  prograde 

•648.63 

044 

-0  41 

182.6  retrograde 

-24.61 

-0.09 

-0  07 

13.7  prograde 

-94  43 

-0.35 

0.01 

13.7  retrograde 

-3.64 

-0.02 

0  04 

365.3  prograde 

25  73 

0.07 

0  15 

365.3  retrograde 

-33.12 

-2  06 

0.33 

27.6  prograde 

14  52 

0.00 

006 

27.6  retrograde 

-13.78 

0.04 

0  00 

121.7  prograde 

-21.58 

-0.10 

-0.03 

121.7  retrograde 

-0.96 

-0  05 

-0  04 

9.1  prograde 

-12.52 

-0.09 

0.00 

9.1  retrograde 

-0.45 

0.01 

-0.02 

31.8  prograde 

3  17 

-0  03 

-0  03 

31.8  retrograde 

-3.02 

0.07 

-0.01 

433.2  prograde 

0.04* 

-0  03 

-0.02 

433.2  retrograde  FCN 

0.33T 

0.32 

-0  07 

*The  standard  deviation  of  each  correction  is  estimated  to  be  0.10  mas,  except 
for  those  of  the  annual  and  433-day  terms  whose  standard  deviations  are  estimated  to 
be  0.12  mas  The  standard  deviations  are  obtained  using  the  techniques  described  in 
Herring  et  al.  ( 1986),  with  the  root-sum-squares  addition  of  a  contribution  of  0  07  mas 
to  account  for  (»)  the  effects  of  the  truncation  of  the  1AU  i960  nutation  series  to  the 
nearest  0  1  mas,  (ii)  the  possible  effects  of  atmospheric  excitation  of  the  nutation,  and 
(iii)  the  effects  of  ocean  and  earth  tide  modeling  errors  in  the  analysis  of  the  VLBI  data 
See  also  Herring  et  al.  (1986)  for  definitions  of  the  quantities  listed. 

tThe  amplitudes  of  these  terms  are  the  root-sum-squares  of  the  amplitudes  of  the 
real  and  imaginary  components. 


that  the  complex  amplitude  of  the  FCN  would  have  a  variance  >(0.5  mas/year)3. 
(Although  formally  this  estimate  of  the  variance  for  the  change  in  the  complex 
amplitude  of  the  FCN  per  year  is  a  lower  bound  (due  to  our  use  of  bounds  in 
its  calculation),  our  model  probably  predicts  too  much  variation  because  of  the 
assumed  white-noise  excitation  which  has  a  constant  power-density,  independent 
of  frequency.)  The  VLBI  estimates  of  the  changes  in  both  the  FCN  and  the 
prograde  signal  are  shown  in  Figure  2.  The  particular  paths  followed  by  each  of 
these  signals  are  not  important.  The  important  point  is  the  absence  of  evidence 
that  the  single  ("average’')  estimate  of  the  amplitude  of  the  FCN  obtained  from 
the  VLBI  data  u  greatly  diminished  by  fluctuation  in  the  complex  amplitude  of 
the  FCN  during  the  six-year  interval  spanned  by  the  VLBI  data.  We  also  note 
that  the  average  value  of  the  amplitude  of  the  FCN  is  consistent  with  our  previous 
•ingle  estimate,  and  much  smaller  than  the  bound  implied  from  the  atmospheric- 
pressure  data. 

The  results  given  in  Table  1  and  Figure  2  indicate  that  the  amplitude  of  the 
FCN  is  no  larger  than  0.6  mas  (99.5%  confidence  limit).  The  VLBI  upper  bound 
for  the  FCN  amplitude  thus  suggests  that  the  PSD  of  the  excitation  is  no  more 
than  0.06  (gm  cm-3)3  day  (99.9%  confidence  limit),  about  five  times  smaller  than 
that  calculated  by  Eubonki  et  al.  Since  we  have  used  an  upper  bound  for  the  FCN 
amplitude  and  a  lower  bound  for  the  damping,  the  actual  PSD  of  the  excitation 
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Figure  1.  Estimates  of  SAi  and  (5A^sinf6  from  the  VLBI  data  analyzed  in  this 
paper.  The  solid  line  in  each  figure  is  computed  from  the  corrections  given  in 
Table  1  and  from  long-period  terms  (periods  >  18.6  years)  which  are  not  given 
in  the  table  because  tneir  values  are  highly  correlated  and  the  estimates  are  not 
considered  reliable. 


might  be  considerably  less  than  the  bound  given. 
4.  CONCLUSION 


Our  analysis  has  placed  an  upper  bound  of  0.6  mas  (99.5  %  confidence  limit) 
on  the  amplitude  of  the  FCN.  If  the  FCN  is  excited  by  surface  loading  with  periods 
near  one  day,  then  this  upper  bound  for  the  FCN  amplitude  combined  with  a 
lower  bound  for  the  damping  time,  allows  an  upper  bound  to  be  placed  on  the 
power-spectral  density  (PSD)  of  the  surface  loading.  This  upper  bound  for  the 
surface  loading  is  0.06  (gm  cm"3)2  day  (99.9  %  confidence  limit)  which  is  about 
five  times  smaller  than  the  estimate  by  Eubanks  et  al.  of  the  PSD  of  the  surface 
loading  due  to  atmospheric  pressure  variations.  These  results  suggest  that  this 
PSD  is  dominated  by  noise  or  that  the  excitation  of  the  FCN  by  surface  loading 
is  not  as  efficient  as  currently  believed  (Sasao  and  Wohr,  1981). 
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Figure  2.  The  estimates  of  the  temporal  changes  in  the  complex  amplitude  of 
the  FCN  and  the  corresponding  prograde  terms,  computed  using  the  Kalmar, 
filter  described  in  the  text.  The  open  squares  show  the  values  on  January  1  for 
years  between  1981  and  1986,  as  well  as  for  the  initial  and  final  epochs  for  the  data 
span.  The  error  bars,  shown  only  at  the  beginning  and  the  end  of  the  time  interval 
spanned  by  the  data,  are  one  standard  deviation,  computed  using  the  algorithm 
discussed  in  the  caption  of  Table  1. 
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DISCUSSION 

Schuh.  Why  have  you  cot  entered  the  22  nutation  coefficient*  u  aoNe  for  parameter*  in  a  VLB! 

•DluUOfi'’ 

Raply  by  Hairing:  We  haven't  implemented  the  •oftw a/e  yet  to  do  that  We  consider  the  nutation  *_n g  1  e 
approach  to  be  more  flexible,  a*  evidenced  by  our  current  paper 


